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Preface 

The  purpose  of  this  study  was  to  determine  the  significance  of 
deviations  from  linearity  for  finite-amplitude  signals  propagating  in 
fluid  tramsmlssion  lines.  The  application  of  this  model  is  to  help 
determine  when  linear  models  of  fluid  transmission-line  transient 
behavior  are  no  longer  applicable. 

The  potential  use  of  this  model  as  an  aid  in  the  design  of  complex 
pneumatic  systems  requires  that  it  run  on  the  moderately  powerful 
computers  readily  available  to  the  design  engineer.  To  this  end,  the 
finite-difference  model  was  constructed  to  keep  computer  resource 
requirements  to  a  minimum.  As  a  design  aid,  the  finite-difference 
model  is  intended  to  supply  data  l>y  which  a  much  simpler,  modal  model, 
in  state-space  format,  can  be  tuned  to  reflect  the  importamt  features 
of  nonlinear  behavior  predicted  by  the  finite-difference  model.  The 
modal  model,  tuned  to  reflect  the  finite-difference  results,  would 
allow  accounting  of  the  nonlinear  effects  in  the  design  of  complex 
fluid  systems  containing  long  transmission  lines. 

Completion  of  this  study  would  not  have  been  possible  without  the 
understanding  and  assistance  of  mv'  doctoral  committee  chairman.  Dr  M. 

E.  Franke.  I  am  also  indebted  to  my  friend,  Kim  Turney,  who  provided 
much  of  the  support  that  helped  me  complete  this  work. 

Wayne  P.  Chepren 
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ABSTRACT 


The  transient  response  of  a  straight,  cylindrical,  fluid 
transmission  line  to  a  finite-amplitude  pressure  input  is  studied 
numerically  to  identify  deviations  from  linear  behavior.  A  nonlinear, 
finite-difference  method  is  developed  to  predict  the  changes  in  the 
shape  of  a  pressure  wave  as  it  propagates.  The  finite-difference 
algorithm  is  also  applied  to  a  linearized  version  of  the  governing 
equations  to  generate  a  basis  against  which  to  identify  nonlinear 
behavior.  The  finite-difference  method  accounts  for  viscous  losses  and 
heat  transfer.  An  algebraic  model  for  calculating  eddy-viscosity  in 
transient  turbulent  flows  is  included. 

Numerical  simulations  of  a  semi-infinite  transmission  line  are 
performed  using  trapezoidal  and  half-sinusoidal  pulses  as  well  as  a 
termlnated-ramp  as  input  pressure  signal  waveforms.  Significant 
differences  between  linear  and  nonlinear  results  are  found,  even  for 
input  signals  with  maximum  pressures  less  than  two  percent  of  the 
initial,  undisturbed  pressure. 

Results  are  also  obtained  for  lines  of  finite  length  with  a 
blocked,  volume  terminated,  or  open  end,  opposite  the  end  where  the 
input  signal  is  applied.  The  trauisient  response  to  a  terminated-ramp 
input  pressure  waveform  is  studied  for  this  case.  Good  agreement  with 
published  results  is  obtained. 


xviii 


An  engineering  model  is  developed  using  the  method  of  modal 
superposition.  The  modal  model  is  modified  to  include 
frequency-dependent  parameters  for  viscous  aund  heat  transfer  losses. 
Results  produced  by  the  modal  model  for  finite-length  lines  agree  well 
with  results  from  the  finite-difference  model  and  published  data.  All 
features  present  in  the  transient  pressure  and  velocity  waveforms 
produced  at  each  end  of  the  simulated  transmission  line  by  the 
finite-difference  model  are  closely  tracked  by  the  modal  model 
results. 


xix 


FLUID  TRANSMISSION  LINE  DYNAMICS  FOR  SIGNALS  OF  FINITE  AMPLITUDE 


I.  Introduction 


Transmission  lines  are  often  present  in  hydraulic  or  pneumatic 
control  and  pressure  measurement  systems.  These  lines  affect  the 
dynamics  of  system  components  by  imposing  loads  on  components  the  lines 
connect.  The  effect  of  the  imposed  loads  must  be  considered  in  order 
to  predict  the  overall  system  response  to  applied  signals. 

Simple  analyses  of  these  systems  allow  transmission  lines  to  be 
treated  as  lumped-parameter  components;  however,  a  line  must  be 
considered  as  a  distributed-parameter  component  to  accurately  predict 
responses  to  applied  signals  containing  frequency  components  with 
wavelengths  significantly  smaller  than  the  length  of  the  transmission 
line.  The  properties  of  the  trauismission  line  are  then  considered  to 
act  on  the  signal  continuously  as  it  travels  through  the  line. 

The  design  of  fluid  control  systems  may  require  accurate  portrayal 
of  nonlinear  effects  in  long  transmission  lines  to  determine  the  input 
signal  shape  needed  to  produce  a  desired  control  action.  Accurate 
fluid  control  in  the  presence  of  signals  of  finite  amplitude  requires 
an  understanding  of  how  nonlinear  effects  act  to  alter  the  signal 
waveform.  The  behavior  of  a  fluid  transmission  line  treated  as  a 
distributed-parameter  component  can  be  described  by  the  Navier-Stokes 
equations.  Numerical  solution  of  these  equations  will  accurately 
predict  transients  produced  by  signals  of  finite  amplitude.  However, 
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numerical  solution  is  expensive  with  respect  to  computer  resources, 
requiring  several  hours  on  a  supercomputer  performing  millions  of 
operations  per  second  to  produce  a  solution.  Designers  of  fluid 
systems  do  not  always  have  ready  access  to  these  computers  auid  must 
rely  on  much  slower  machines.  A  finite-difference  model  requiring  only 
modest  computer  resources  and  based  on  the  simplified  Navier-Stokes 
equations  would  be  useful  to  fluid  system  designers. 

Despite  the  availability  of  a  finite-difference  model  capable  of 
running  on  computers  easily  accessible  to  designers,  a  simplified 
lumped  parameter  model  capable  of  interfacing  with  standard  techniques 
used  in  system  design  is  also  needed.  Numerical  solution  of  the 
partial  differential  equations  governing  fluid  flow  results  in  an 
overwhelming  amount  of  data  that  can  be  cumbersome  and  expensive  to 
apply  to  complex  fluid  systems.  Lumped-parameter  models  have  been  used 
for  many  years  to  calculate  the  transient  response  of  structures  to 
suddenly  applied  or  nonperiodic  loads.  These  techniques  are  based  on 
modal  superposition  which  is  mathematically  valid  for  linear  systems 
only.  However,  modifications  to  modal  superposition  that  account  for 
nonlinear  effects  in  structures  have  produced  good  results  in  some 
cases  (47:72).  Additionally,  state-space  analyses  of  physical  systems 
have  relied  on  describing  functions  to  modify  linear  models  to  account 
for  nonlinear  effects  in  transient  solutions. 

Background 


Accurate  predictions  of  the  transient  behavior  of  the  fluid 
systems  have  been  obtained  using  modal  approximations  to  invert  the 
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we 11 -developed  small  laminar  disturbance  (SLD)  models  from  the  Laplace 
domain  to  the  time  domain  (24).  Due  to  the  linear  nature  of  the  SLD 
models,  the  modal  approximations  cannot  be  modified  to  deal  with  a 
turbulent  mean  flow,  signals  of  finite  amplitude,  or  transmission  line 
terminations  with  complicated  nonlinear  behavior.  Other  transient 
solutions  based  on  a  SLD  model  suffer  from  similar  drawbacks  (5:78;  7; 
30). 

A  simplified  engineering  model  capable  of  predicting  the  effects 
of  f Inite-aunplitude  signals  and  viscous  losses  for  transient  turbulent 
flow  in  fluid  transmission  lines  should  be  useful  in  pneumatic  system 
design.  This  engineering  model  should  be  compatible  with  staindard, 
control -system  analysis  and  design  techniques.  It  should  also  require 
only  a  minimum  of  empirical  Information  to  model  transmission- line 
response  and  be  compatible  with  ordinary  differential  equations  used  to 
describe  other  system  components.  Such  a  model  could  be  developed  from 
solutions  generated  with  a  finite-difference  model.  The  engineering 
model  will  then  contain  parameters  whose  values  are  adjusted  to  give 
its  solutions  the  same  importamt  characteristics  as  solutions  produced 
by  the  computer  intensive  finite-difference  model.  The  adjustable 
parameters  of  the  engineering  model  allow  the  information  obtained  with 
the  f *nite-difference  model  to  be  compressed  into  a  convenient  format. 

Ideally,  results  of  the  numerical  model  obtained  for  a  small 
number  of  cases  could  be  used  to  modify  the  engineering  model  for  a 
wide  range  of  signals  and  system  configurations.  Although  the  need  for 
actual  system  construction  and  testing  will  never  be  eliminated,  the 
finite-difference  and  engineering  models  used  together  in  the  early 
stages  of  design  should  result  in  a  prototype  fluid  system  with 
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improved  performance.  Less  redesign  and  rework  of  the  fluid  system 
during  protot'T>e  testing  would  be  required  to  achieve  desired 
performance. 

Previous  Investigations 


Propagation  of  sinusoidal  signals  through  fluid  transmission  lines 
with  no  established  mean  flow  has  been  studied  by  numerous 
investigators  using  frequency-domain  techniques  (26;  29;  49;  55).  The 
assumption  of  a  SLD  allows  the  continuity,  momentum,  energy,  and  state 
equations  to  be  greatly  simplified  and  solved  to  produce  explicit 
expressions  for  distributed  friction  and  heat  transfer  losses  as  a 
function  of  frequency.  These  frequency-domain  models  can  be  inverted 
from  the  Laplace  domain  to  the  time  domain  to  obtain  the  transient 
response  to  the  simple  cases  of  impulse  and  step  inputs  (7).  Methods 
to  obtain  the  trainslent  response  for  a  suddenly  applied  sinusoidal 
signal  (28'  as  well  as  the  approximate  response  for  an  arbitrary  signal 
(30)  were  developed  assuming  a  small  laminar  disturbance. 

Two  other  approaches  that  rely  on  frequency-domain  models  to 
obtain  the  transient  response  to  an  arbitrary  signal  are  the  frequency 
response  conversion  method  (33;  76;  77;  78;  79;  80;  81;  and  the  Fast 
Fourier  Transform  (FFT)  method  (71).  Frequency  response  conversion  is 
the  direct  conversion  of  the  exact  frequency  response  model  for 
terminated  transmission  lines  into  the  transient  response  (33:3).  The 
fast  Fourier  transform  method  calculates  the  effect  of  the  transmission 
line  on  each  frequency  component  of  the  transformed  signal  and  then 
reconstructs  the  signal  from  the  modified  frequency  components. 
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Methods  to  incorporate  the  effects  of  laminar  mean  flow  into  SLD 
models  in  the  frequency  domain  have  been  developed  by  various  authors 
(5;  51;  52).  In  this  case  the  signal  is  imposed  on  eui  otherwise  steady 
laminar  flow.  The  more  complex  problem  of  accounting  for  the  effects 
of  turbulent  mean  flow  in  the  SLD  models  has  also  been  addressed  19; 

19;  34;  43;  44).  Viscous  losses  in  transient  turbulent  flow  have  also 
been  modeled  (74).  In  addition,  approximations  to  the  transmission 
line  transfer  function  as  well  as  the  complicated  impedance  and 
admittemce  expressions  have  been  made  (25;  50)  allowing  the  transient 
response  for  a  variety  of  input  signals  to  be  easily  calculated. 

With  the  exception  of  the  FFT  method  (71)  and  the  block  diagram 
reduction  method  (30)  the  transient  response  for  arbitrary  signals 
cannot  be  obtained  from  the  SLD  models.  Additionally,  the  use  of  SLD 
models  recpiires  line  terminations  to  behave  in  a  linear  fashion. 
Numerical  techniques  such  as  the  method  of  characteristics  and 
finite-difference  methods  can  be  applied  to  overcome  these  limitations. 

The  method  of  characteristics  is  a  numerical  method  long  used  for 
the  solution  of  transient  flow  in  pipelines  (1;  64:20).  This  method 
has  also  been  applied  to  systems  with  finite-amplitude  signals  in 
compressible  fluids  (3;  41;  60).  A  one-dimensional  finite-difference 
solution  retaining  the  convective  terms  in  the  governing  equations  has 
also  been  presented  (19).  In  this  method,  and  the  method  of 
characteristics,  the  viscous  losses  are  calculated  using  the 
Darcy-Weisbach  friction  factor  developed  for-  steady  flow.  It  has  been 
shown  that  use  of  friction  factors  developer  for  steady  flow 
significantly  underestimates  the  attenuation  of  a  propagating  signal 
(23:6;  74:866).  The  method  of  characteristics  has  been  adapted  to 


handle  viscous  losses  in  transient  flow  by  developing  forcing  functions 
based  on  earlier  frequency  domain  friction  models  (9;  82).  This  method 
works  well  for  lines  with  laminar  mean  flow  and  for  a  limited  range  of 
turbulent  mean  flows  (69). 

Solutions  of  the  full  Navier-Stokes  equations  and  of  the 
axisymmetric  equations  for  pulsed,  turbulent  pipe  flow  have  been 
obtained  (35;  54).  These  numerical  approaches  provide  detailed 
information  on  the  transient  behavior  of  the  flow,  however  a  large 
number  of  grid  points  in  the  radial  direction  is  required  to  accurately 
represent  the  boundary- layer  phenomena  (70:588). 

A  significant  reduction  in  the  number  of  grid  points  may  be 
realized  by  simplifying  the  Navier-Stokes  equations  to  produce  a  set  of 
equations  that  treat  heat  transfer  and  friction  terms  as  forcing 
functions  (70).  By  using  piecewise  defined  polynomials  to  approximate 
velocity  and  temperature  profiles,  the  number  of  mesh  points  needed  to 
approximate  the  friction  and  heat  transfer  forcing  functions  is 
significantly  decreased.  The  availability  of  the  velocity  and 
temperature  profiles  from  the  polynomials  allows  the  use  of 
Boussinesq-type  friction  models  to  calculate  the  values  of  the  damping 
functions. 

A  numerical  method  that  uses  modal  superposition  to  predict  the 
transient  response  of  a  transmission  line  has  also  been  presented  (31; 
40).  Modal  superposition  removes  spatial  dependence,  producing  a 
system  of  ordinary  differential  equations  compatible  with  standard 
methods  of  system  ainalysis. 

The  studies  discussed  above  are  categorized  in  Table  1  according 
to  the  solution  method  provided  by  the  model.  A  more  detailed  summary 
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of  the  frequency-domain  models  has  been  presented  by  Goodson  and 
Leonard  (22). 

Objective 


The  intent  of  this  study  is  to  determine  the  importance  of 
nonlinear  effects  in  the  transient  response  of  a  fluid  transmission 
line  for  signals  of  finite  amplitude.  An  explicit,  nonconservative, 
finite-difference  model  accounting  for  nonlinear  effects  as  well  as 
viscous  and  heat  transfer  losses  in  tramsient  flow  is  developed.  The 
model  has  provisions  for  both  laminar  and  turbulent  viscous  losses  and 
heat  transfer  in  addition  to  provisions  for  a  mean  flow.  The  transient 
response  of  blocked,  open,  and  volume-terminated  transmission  lines  is 
examined  to  Identify  nonlinear  effects  and  assess  limitations  of  the 
finite-difference  model  for  this  application.  A  linear  version  of  this 
model  is  also  developed  to  compare  linear  and  nonlinear  transmission 
line  behavior.  The  relative  importance  of  viscous  and  heat  transfer 
losses  is  also  examined.  The  finite-difference  model  is  developed  to 
minimize  the  computational  capability  required.  The  intent  is  to  allow 
the  finite-difference  model  developed  in  the  present  study  to  run  on 
computers  of  modest  capability  and  find  use  as  a  fluid  system  design 
tool . 

An  engineering  model  suitable  for  use  with  the  state-space  method 
of  dynamic  system  analysis  is  also  developed.  This  model  allows  the 
influence  of  the  transmission  line  to  be  accounted  for  in  the  transient 
response  for  a  complete  fluid  system. 
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TABLE  1 


Transmission  Line  Literature  Summary 


Invarslon  of  fraquancy  donaln 

Fraquancy  doaa  In  aolutlona 

Bodals  to  got  translont 

ao 1 ut Ion 

Ibarall  (19S0) 

Brown  (1962) 

Rolaaann  and  Grogan  (19S7) 

Kantola  (1971) 

Nlchola  (1961) 

Karan  (1972) 

Broun  (1962) 

Karan  and  Laonard  (1973) 

Karan  and  Franka  (1967) 

Fraquancy  raaponaa 
convarslon 

Katz  (1976) 

Yahloka  (1983) 


Faat  Fourier 

tranafom  nathod 

Tsang  at  al 

(1984/ 

Fraquancy  dona  In  frlctlo- 
nodals  Incorporating 
affacta  of  thrrughflow 

Approxlnatlona  to  tranafar 
natrlx  and  1  npadence  tarna 

Brown  at  al  (1969) 

Oldanburgar  (1964) 

0.nar(1969) 

Hul lender,  Hooda, 

Hood  and  Funk  (1970) 

Narcler  and  Hrlght  (1973) 
Katz  at  al  (1974) 

Bradan  (1973) 

Margolla  and  Brown  (1975) 

Platclier  (1979) 

Brown  (1981) 

and  Hau  (1983) 

Mathod  of  charac t rjr  1  s 1 1  cs 


Banson  at  al  (1964) 
Manning  (1968) 

Sandoz  (1973) 


Method 

of  characteristics 

with  unsteady  friction 

Z 1  e 1 ke 

(1968) 

Ntasar leal  app roaches 

One  dimensional  nunerlcal 

(Two  or  three  dloanslonal) 

method 

Tsao  (1968) 

Funk  and  Robe  (1970) 

Klrmse  (1979) 

Reddy  at  al  (1985) 
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Approach 


The  finite-difference  model  will  be  based  on  the  approach 
described  by  Tsao  (70).  The  governing  equations  are  simplified  and 
expressed  in  what  Tsao  calls  damped-wave-pi  opagation  equations 
(70:589).  These  damped-wave  equations,  as  they  will  be  called  in  this 
study,  will  then  be  solved  using  the  Gabutti  finite-difference 
algor itnm  (<d:294;  21)  with  an  algebraic  turbulence  moael  for  viscoua 
losses  in  trcinsient  turbulent  flow  that  approximates  the  apparent 
viscosity  at  turbulent  Reynolds  numbers.  An  algorithm  for  generating 
the  pressure  and  velocity  fields  for  an  established  mean  flow  for  use 
with  the  finite-difference  model  is  also  developed.  The  modal 
superposition  method  described  by  Karnopp  (31)  and  Lebrun  (40)  is  used 
to  develop  a  linear  engineering  model. 

Organization  of  Thesis 

The  complete  Navler-Stokes  equations  are  reduced  to  a  simplified 
form  in  Chapter  II  using  the  well-known  transmission  line  assumptions. 
Following  Tsao,  the  simplified  equations  and  an  equation  of  state  are 
combined  to  produce  a  system  of  nonlinear  equations,  called  damped-wave 
equations  (70).  The  use  of  piecewise  defined  polynomial  approximations 
of  the  velocity  zuid  temperature  profiles  is  introduced  and  applied  to 
the  solution  of  the  damped-wave  equations. 

Damping  functions  that  are  used  with  the  damped-wave  equations  and 
the  piecewise  polynomial  approximations  developed  in  Chapter  II  are 
presented  in  Chapter  III.  Development  of  a  relatively  simple, 
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algebraic  turbulence  model  suitable  for  transient  flows  and  a  method 
which  calculates  the  initial  pressure  and  velocity  field  for  a 
turbulent  or  laminar  mean  flow  are  also  shown  in  Chapter  III.  The 
Gabutti  numerical  algorithm  used  to  solve  the  damped-wave  equations  is 
discussed  in  Chapter  IV.  Features  of  other  numerical  methods  are  also 
discussed  and  compared  to  those  of  the  Gabutti  algorithm. 

An  engineering  model  based  on  modal  superposition  auid 
electromagnet !r  transmission  line  theory  is  developed  to  apprsxiactc 
the  finite-difference  solution  of  the  damped-wave  equations  in 
Chapter  V.  A  method  for  determining  the  modal  parameters  from  linear 
transmission  line  theory  is  described. 

Numerical  results  from  both  models  are  presented  amd  discussed  in 
Chapter  VI.  Results  are  also  compared  with  results  from  past  studies 
where  possible. 

Chapter  VII  contains  conclusions  and  recommendations  for  further 
studies. 
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I I .  Damped  Wave  Equations 


A  finite-difference  model  capable  of  accurately  representing  the 
dynamic  response  of  a  fluid  transmission  line  to  flow  and  pressure 
disturbances  applied  at  the  ends  of  the  line  is  developed  in  this 
study.  This  numerical  method  must  accoiont  for  nonlinear  effects  due  to 
finite-amplitude  signals,  effects  due  to  the  presence  of  an  initial 
turbulent  or  laminar  mean  flow,  and  both  laminar  or  turbulent  viscous 
losses  and  heat  transfer.  The  development  of  this  finite-difference 
model  can  be  separated  into  three  steps;  simplification  of  the 
aixisymmetric  Navier-Stokes  equations  into  a  manageable  form, 
development  of  suitable  expressions  for  viscous  losses  and  heat 
transfer,  and  selection  and  application  of  the  numerical  algorithm  to 
be  used  to  solve  the  system  of  equations  resulting  from  the  first  two 
steps.  This  chapter  describes  the  first  step  in  development  of  the 
finite-difference  model. 

The  simple  geometry  of  a  cylindrical  transmission  line  allows 
assumptions  that  greatly  simplify  the  Navier-Stokes  equations. 

Following  common  practice,  the  terra  Navier-Stokes  will  include  the 
continuity  and  energy  equations  as  well  as  the  components  of  the 
viscous  momentum  equations. 

The  following  assumptions  can  be  made  for  the  transmission  line 
shown  in  Figure  1. 

1.  The  transmission  line  is  rigid,  circular  in  cross  section,  and  has 
constant  cross-sectional  area. 

2.  The  gas  behaves  ideally.  Liquids  are  assumed  to  be  barotropic 
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Fig.  1.  Partitioning  of  Flow  into  Four  Intervals 


(p  =  p(p))  with  fluid  properties  such  as  the  speed  of  sound  and 
the  ratio  of  specific  heats  assumed  continuous. 

3.  The  mean  flow  is  fully  developed  leuninar  or  turbulent  flow. 

4.  The  time  varying  pressure  is  uniform  across  any  given  cross  section 
of  the  line. 

5.  The  flow  is  at  all  times  parallel  to  the  axis  of  the  line.  This 
means  that  velocity  components  in  the  0  and  r  directions  have  zero 
magnitude. 

6.  The  transmission  line  walls  are  isothermal. 

7.  The  length  of  the  transmission  line  is  much  greater  thaui  the 
radius.  This  allows  the  end  effects  and  higher  mode  propagation 
to  be  neglected. 

8.  The  bulk  viscosity  is  assumed  to  be  zero  and  the  fluid  undergoes  no 
relaxation. 

9.  There  are  no  significant  body  forces  acting  on  the  fluid. 

A  set  of  nonlinear,  damped-wave  equations  may  be  obtained 

following  a  procedure  similar  to  that  used  by  Tsao  (70).  The 
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nonlinear,  damped-wave  equations  with  the  convective  terms  retained  are 
shown  below.  Details  of  the  development  of  these  equations  are  shown 
in  Appendix  B: 

du  c^a(Ln[p])  du  c^  1  3  — 

—  + - +  u —  = - [r(T  +  puv)]  (1) 

at  7  dx  dx  p7  r  dr 

a(Ln[pl)  au  a(Ln(pl)  (y-l)l  a 

-  +  7 —  +  u -  - - [r(q  +  pvc  T  ]  (2) 

at  dx  dx  p  r  dr 

The  right-hand  side  of  Eq  (1)  represents  the  viscous  losses  while  the 
corresponding  term  in  Eq  (2)  represents  losses  due  to  heat  transfer. 

All  variables  in  Eqs  (1)  and  (2)  are  time  averaged  except  where  noted 
by  primes.  Primes  denote  values  that  fluctuate  about  the  average.  The 
set  of  initial  and  boundary  conditions  are 


p  =  f(x),  u  =  g(x,r),  T-T^  =  0 

(p)  =  f(t),  or  (u)  =  f(r,t) 

x=0  x=0 

(u)  =  g(r, t)  or  (p)  =  f (t) 

x=L  x=L 


'  au  ' 

'  aT  ' 

.  3r  ^ 

r  =  0 

.  Sr  . 

(u)  =  (T-T  )  =  0 

r=a  H  r=a 


t<0,  OsxaL,  Osrsa 
Ost,  Osrsa 
Ost,  Osrsa 

Ost,  OsxsL 

0  £  t,  OsxsL 


The  above  equations  are  nondiraensionalized  using  the  groupings 
shown  below: 


u  = 


P  = 


2 

p  C 
0  0 


X  = 


L 


L 


T  = 


T 


u 


r 

a 
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Before  doing  any  further  work  with  the  damped-wave  equations,  it 
is  necessary  to  develop  a  means  of  specifying  the  right-hand  side  of 
Eqs  (1)  and  (2).  These  expressions  will  be  referred  to  as  damping 
functions,  since  they  represent  attenuation  resulting  from  viscous  and 
heat  transfer  losses.  The  method  proposed  by  Tsao  uses  averaged 
velocity  and  temperature  functions  as  a  means  of  approximating  these 
damping  functions  (70:591).  A  conventional  differencing  scheme 
requires  a  very  fine  mesh  near  the  pipe  wall  in  the  radial  direction  to 
represent  accurately  the  behavior  of  the  velocity  profile  in  the 
bovindary  layer.  The  averaged  velocity  and  temperature  functions  in 
conjunction  with  piecewise  defined,  polynomials  are  used  to  generate 
approximate  velocity  and  temperature  profiles  over  the  cross  section, 
significantly  reducing  the  number  of  mesh  points  required  in  the  radial 
direction. 

The  cross  section  of  the  transmission  line  is  partitioned  into  M 
flow  intervals  as  shown  in  Figure  1.  The  locations  of  the  interval 
boundaries  with  respect  to  the  centerline  are  denoted  as 
0=r  <r  <r  <***<r  =a  where  a  is  the  internal  radius  of  the 

0  12  M 

transmission  line.  The  area  of  the  annular  region  formed  between  r 

n*l 

2  2 

and  r  is  denoted  by  AA  and  is  defined  as  AA  =  ti  (r  -  r  )  where 

n  n  n  n  n-1 

n  =  1,2 . M.  The  averaged  temperature  and  averaged  velocity 

functions  are  defined  as 

r 

1  r " 

T  (x,t)  =  -  T(x,  r ,  t )2Trrdr  (3) 

"  A  A  J 

AA  r 
n  n-1 
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If" 

u  (x, t)  =  -  u(x, r, 

"  AA  J- 


t )  Zrrrcir 


(4) 


AA  -r 
n  n-1 


These  averaged  functions  are  then  Incorporated  into  the  damped-wave 
equations  to  produce  a  set  of  simultaneous,  one-dimensional,  partial 
differential  equations: 


^  c^a(Ln[p]) 

au 

+  U-"  = 

at  7  ax 

"ax 

2irc^  r 
ypAA  L 

n 

'’n 

n 

a(Ln[p]) 

^  ..  a(Ln(p]) 

at 

-  ax 

n  =  1,  2, 


Ju 


ax 


a  p 


fail 


(5) 


.  M 


(6) 


where  u  =  - f  uZitrdr  and  AA  =  f  Znrdr  ,  auid  p  and  k  are  the 

.  J  J  T  T 

AA  ''o  •'o 

total  fluid  viscosity  and  thermal  conductivity,  respectively.  Total 
fluid  viscosity  is  understood  to  be  the  sura  of  the  molecular  fluid 
properties  and  the  apparent  fluid  properties  due  to  turbulence,  and 
are  found  using  a  suitable  turbulent  viscosity  model  that  will  be 
described  in  the  next  chapter.  The  groupings  shown  below  are  used  to 
nondimenslonalize  p  and  k  : 

T  T 


k  =  k 

T  T 


r  T  L  1 

w 

• 

f  1  S  If 

o 

o 

r 

2 

c  p  a 

0  °  J 
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The  energy  equation,  in  terms  of  averaged  variables,  is  written  as 


31 

I 

at 


- 

+  u  — " - — (Ln[p] )  - 

"ax  Tc  at 

p 


2c 


,2  2  . 

yc  p  (r  -  r  ) 

p  n  n-1 


r  k  (r  ) 

n  T  n 


a(LnIp] ) 
yc  "  ax 

P 


'  ar  ' 

-  r  k  (r  ) 

ar 

V.  > 

n-1  T  n-1 

r=r 

n 

ar 

k.  J 

r=r 

n-l-J 

n  =  1,  2, 


(7) 

••  M 


The  terms  — ^  and  u  ,  required  to  numerically  solve  the  energy 
equation,  are  obtained  from  solution  of  the  damped-wave  equations. 

The  derivatives  of  velocity  and  temperature  at  the  interval  boundaries 
are  necessary  to  calculate  values  for  the  right-hand  side  of  Eqs  (5) 
through  (7).  Tsao  describes  a  method  for  reconstructing  these  profiles 
(70:595).  The  reconstruction  uses  piecewise  f>olynomlals  and  the  values 
of  the  averaged  variables  u  auid  T  in  a  mamner  similar  to  a  cubic 

n  n 

spline  function  to  approximate  the  velocity  amd  teu..  jrature  profiles. 
Let  each  interval  be  denoted  by  I  with  n  =  1 ,2, • • • ,M  as  shown  in 

n 

Figure  1.  In  each  interval,  the  velocity  is  approximated  by  a 
polynomial  (70:591): 


u(x,r,  t)3-{ 


P  ( 

n 


2 

x,r,t)=a+br+cr  rel 

111  1 

2  3 

x,r,t)=a  +b  (r-r  )+c  (r-r  )  +  d  (r-r  ) 

n  n  n-1  n  n-1  n  n-1 

rel,  n=l,  2, 


(8) 


Equation  (8)  has  4M-1  unknown  coefficients  which  require  4M-1 
linearly  Independent  equations  to  determine.  The  value  for  the 
averaged  velocity  in  each  interval  is  known  from  the  initial  conditions 
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or  from  the  solution  of  Eq  (5)  at  each  time  step.  The  definition  of 
averaged  velocity  shovm  in  Eq  (5)  is  applied  to  Eq  (8)  and  used  with 


the  known  values  of  averaged  velocity  to  provide  M  independent 
equations.  Each  piecewise  polynomial  as  well  as  its  first  and  second 
derivatives  are  required  to  be  continuous  at  each  boundary  between 
intervals.  This  provides  an  additional  3M-3  independent  equations. 

Two  more  equations  result  from  the  assumptions  of  no  slip  at  the  wall 
and  axisymmetric  flow  where  u)  and  are  zero,  respectively. 

r=a  or  J  r  =  0 

Each  of  the  4M-1  unknowns  in  the  polynomials  is  determined  at  each  time 
step  from  a  set  of  recursion  relationships  resulting  from  the 
application  of  these  conditions.  The  procedure  for  determining  these 
recursion  relationships  is  shown  in  detail  in  Appendix  C.  The 
recursion  relations  are  used  to  produce  a  system  of  equations  in  terms 
of  the  coefficients  c  and  d  of  Eq  (8).  Solution  of  this  system  of 

n  H 

equations  allows  reconstruction  of  the  velocity  profiles  and  gradients: 


- 1 


c 

1 

1 

1 

1  1 

u 

2 

1.1 

1.2 

1,M-1 

l.M 

1 

c 

1 

1 

1 

1 

U 

3 

• 

2.1 

• 

2,2 

• 

2,  M-1 

• 

2,M 

• 

2 

• 

• 

= 

• 

•  •  •  • 

• 

• 

• 

0 

C 

1 

0 

1 

• 

1 

0 

1 

• 

U 

N 

M-1, 1 

M-1, 2 

M-1,M-  1 

M-1,M 

M-1 

d 

1 

1 

1 

1 

M 

M,  1 

M,2 

M,M-1 

M,M 

U 

M  -• 

The  elements  1  of  Eq  (9)  are  constants  determined  by  the  locations 

i.  J 

of  the  interval  boundaries  r 

n 

Temperature  profiles  are  determined  in  the  same  manner  except  that 
the  profile  of  T  -  T^  is  found  instead  of  the  actual  temperature 


profile,  with  T  -  T  replacing  u  and  T  -  T  replacing  u  in  Eqs  (8)  and 
(9).  is  the  wall  temperature  of  the  trsuismission  line  and  remains 

constant.  The  field  of  values  of  u  is  provided  by  the  numerical 

n 

solution  of  Eqs  (5)  through  (7)  from  which  Eq  (9)  produces  values  of  c 

n 

and  d  at  each  time  step.  These  values  are  substituted  into  the 
recursion  relationships  for  a  ,  b  ,  amd  d  (found  in  Appendix  C), 

n  n  n 

Q 

allowing  ^^u(JAx,r,t)  to  be  determined  from  Eq  (8).  Values  of  the 
velocity  gradients  at  r  =  r  are  used  to  evaluate  the  right-hauid  side 

n 

of  Eq  (5).  This  entire  sequence  is  repeated  at  each  time  step 

throughout  the  numerical  solution  of  Eqs  (5)  through  (7).  In  exactly 

the  same  manner,  -5—  T(JAx,r,t)  is  evaluated  from  values  of  T  -  T 
or  n  w 

produced  by  the  numerical  solution.  The  right-hand  sides  of  Eqs  (6) 
and  (7)  are  then  evaluated.  This  procedure  must  include  an  additional 
step  when  the  flow  is  turbulent.  Expressions  for  the  total  viscosity 
amd  the  total  thermal  conductivity  must  be  evaluated  at  each  time 
step. 

The  use  of  piecewise  defined  polynomials  to  approximate  the 
velocity  and  temperature  profiles  amd  gradients  is  very  similar  to 
cubic  spline  interpolation.  Cubic  spline  interpolation  provides  a  goed 
approximation  of  slopes  which  has  made  it  popular  as  a  means  of 
numerical  differentiation. 

A  linearized  set  of  equations  obtained  from  the  Navier-Stokes 
equations  was  presented  by  Tsao  (70:589).  Using  these  equations  and 
proceeding  as  described  in  Appendix  B,  the  following  equations  are 
obtained: 
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fH"  +  c^aCLnfpJ) 

3t  y  ax 
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ypAA  L  n  ^ 

^  n  n-l  ^ 

''  n-l-* 

n  =  1,  2,  •••,  M  (10) 


a(Ln[p] )  ^ 
at 


au 
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'af 

ax 

a  p  ^ 

ar 

k  y 
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_”  +  _2_  ^(Ln(p])  = 

at  yc  at 
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2c‘ 


yc  p  (r^-  ) 

p  n  n-l 


r  k  (r  ) 

'  aT  ' 

-  r  k  (r  ) 

f  > 

aT 

n  T  n 

ar 

k.  i 

n“i  1  II*  1 

r=r 

n 

ar 

1  > 

r=r 

n-lJ 

n  =  1,  2,  •••,  M  (12) 


These  equations  are  in  the  s.-'me  format  as  Eqs  (5)  through  (7).  The 
numerical  solutions  of  both  sets  of  equations  will  be  compared  in 
Chapter  VI  to  identify  nonlinear  effects. 
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III.  Damping  Functions 


Model  for  Viscous  Losses  in  Transient  Turbulent  Flow 

Data  on  the  transient  turbulent  flows  under  consideration  in  this 
study  is  very  limited.  As  a  result,  no  criteria  are  available  on  which 
to  base  selection  of  a  turbulent  viscosity  model.  A  simple,  algebraic 
model  to  approximate  of  the  effect  of  viscous  losses  due  to  turbulence 
in  the  solution  of  the  equations  developed  in  Chapter  II  will  be 
discussed  in  this  chapter. 

Systems  of  interest  in  this  study  include  the  presence  of  a  mean 
flow.  This  is  a  steady,  fully-developed  flow  assumed  to  be  in 
existence  before  application  of  any  signal.  Determination  of  the 
initial  velocity  field  for  a  mean  flow  would  normally  require 
continuing  the  numerical  solution  of  the  damped-wave  equations 
developed  in  the  previous  chapter  until  steady  state  is  reached.  An 
algebraic  method  using  the  turbulent  viscosity  model  to  generate  a 
velocity  field  for  use  as  an  initial  condition  will  be  developed  in 
this  chapter. 

Expressions  for  the  total  viscosity  and  the  total  thermal 
conductivity  are  needed  to  solve  Eqs  (5)  through  (7)  and  (10) 
through  (12).  For  the  case  of  laminar  flow,  and  k^  are  identical  to 
fi  and  k  respectively.  The  presence  of  turbulence  complicates  the 
problem  significantly. 

Viscous  effects  resulting  from  turbulent  flow  are  often  treated  in 
a  quasi-steady  manner  when  numerical  techniques  such  as  the  method  of 
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characteristics  are  applied  to  fluid  transmission  line  problems.  In 
this  case,  viscous  forces  are  calculated  using  expressions  valid  only 
for  fully  developed  velocity  profiles.  In  the  quasi-steady 

pU^fL 

approximation,  the  friction  force  is  usually  expressed  as  Ap  = 
where  f  is  the  Darcy-Weisbach  friction  factor  and  is  calculated  as  a 
function  of  the  Reynolds  number.  Examples  of  this  approach  are 
numerous  (3:117;  16:9;  75:18). 

The  effectiveness  of  this  quasi-steady  approach  for  representing 
viscous  effects  resulting  from  turbulent  flow  has  been  questioned.  For 
example.  Wood  and  Funk  found  that,  for  the  extreme  case  of  pressure 
wave  propagation  in  a  tube,  the  quasi-steady  approximation  fails  to 
predict  the  true  attenuation  of  the  wave  (74:866).  They  use 
waterhammer  caused  by  sudden  valvet  closure  as  an  exeunple  and  show  that 
the  large  reverse  pressure  gradient  reverses  the  flow  near  the  wall, 
producing  significant  dissipation.  This  dissipation  decreases  the 
magnitude  of  the  wave  more  than  a  quasi-steady  calculation  would 
indicate. 

Expressions  for  viscous  losses  in  transient  flow  of  liquids  based 
upon  the  frequency  dependent  laminar  friction  models  discussed  in 
Chapter  I  have  been  developed  for  laminar  flows  (82)  amd  modified  for  a 
limited  range  of  turbulent  flows  (69).  These  techniques  use  weighting 
functions  to  calculate  damping  based  on  the  time  rate  of  change  of  the 
average  velocity  and  the  velocity  history.  This  method  works  well  for 
laminar  flows  and  when  steady,  turbulent,  mean  flows  are  present.  The 
methods  of  Zielke  (82)  and  Trikha  (69),  however,  are  one  dimensional 
and  provide  no  means  to  calculate  viscous  losses  from  transient 


turbulent  flow. 


The  equations  developed  In  Chapter  II  produce  two-dimensional 


velocity  and  temperature  profiles  which  allow  the  use  of 
Bousslnesq-type  eddy  viscosity  models  to  predict  the  viscous  and  heat 
transfer  losses  on  disturbances  traveling  through  the  transmission 
line.  A  variety  of  eddy  viscosity  models  commonly  used  In  numerical 
methods  may  be  employed  to  caJ''”lat«  values  for  and 

Launder  defines  a  turbulence  model  as  a  set  of  differential  or 
algebraic  equations  whose  solution  enables  certain  features  of 
turbulent  flow  to  be  predicted.  Classifications  of  turbulence  models, 
as  described  by  Launder  (39),  are  shown  In  Table  2. 


Table  2 

Classification  of  Turbulence  Models  (39:417) 


Type 

Number  of  Dependent 

Turbulent 

of  Model 

Transport  Equs. 

Variables 

Mixing  Length/ 

Eddy  viscosity 

0 

— 

Bradshaw’ s  Method 

(also  Nee  and  Kovasznay) 

K 

Two  equation  Viscosity 

Model  2 

k,  E 

or 

equivalent 

Stress  Transport  Model 

5-7 

/  / 

u  u  ,  c 

1  J 

Third  Order  Closure 

9-17  u  V  , 

1  J 

/  /  / 

U  U  U  ,  G 

1  J  k’ 

Boussinesq  suggested  that  an  apparent  scalar  eddy  viscosity  might 
relate  turbulent  shearing  stresses  to  the  rate  of  mean  stress  (61:189). 
Using  an  analogy  with  the  kinetic  theory  by  which  the  molecular 
viscosity  for  gases  has  been  calculated,  the  turbulent  viscosity  is 
modeled  as  u  =  pv  I  where  v  and  t  are  characteristic 

xurb  Turb  Turb 

velocity  and  length  scales  for  turbulence,  respectively  (2:221).  The 

first  three  models  listed  in  Table  2  provide  a  means  to  solve  for  v 

and  1.  The  remaining  two  models  are  solved  for  the  Reynolds  stresses 

directly  but  require  the  solution  of  several  simultaneous,  differential 

equations.  The  averaged  velocity  functions  developed  in  Chapter  II 

allow  the  use  of  one  of  the  first  three  types  of  turbulence  models 

shown  in  Table  2  to  obtain  the  values  of  eddy  viscosity. 

Eddy  viscosity  models  were  first  developed  for  steady  flow 

applications.  The  goal  of  these  models  is  to  find  a  relationship 
/  / 

^Turb  *  ~  ^(J«»r,u,p)  that  will  faithfully  reproduce  the  solid 

curve  through  the  data  on  the  universal  wall  plot  for  steady,  turbulent 
boundary  layers  on  smooth,  solid  surfaces.  From  an  analogy  with 
laiminar  flow,  the  eddy  viscosity  formulation  of  Boussinesq  has  the  form 
(73:559) 


T 

Turb 


/  / 


puv  =  fl 


Turb 


au 

dy 


Prandtl  suggested  the  mixing  length  formulation  (73:469) 


(13) 


T 

Turb 


-puv  =  pt 


ay 


ay 


(14) 
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where  t  Is  some  effective  interaction  distance  between  eddies.  The 

■ 

eddy  viscosity  and  the  mixing  length  formulations  are  essentially 
equivalent  since  we  may  write  from  Eqs  (13)  and  (14)  as 


H  ~ 
'^Turb 


ay 


(15) 


It  is  important  to  note  that  the  turbulent  viscosity,  .  must  be 

expected  to  depend  on  the  state  of  the  flow  aind  not  Just  the  state  of 
the  fluid. 

The  approach  above  is  known  as  a  gradient  transport  formulation 

and  produces  t  =0  where  the  velocity  profile  has  zero  slope. 

Turb 

Due  to  the  symmetry  of  the  flows  considered  in  this  study,  there  is  no 
net  momentum  transfer  across  the  centerline  and  the  gradient  transport 
formulation  applies.  Schetz  points  out  that  experiments  with 
high-velocity  injection  into  a  turbulent  boundary  layer  through  an 
upstream  slot  in  the  wall  indicate  that  a  nonzero,  turbulent  shear 
stress  may  exist  at  local  maoclma  or  minima  in  the  velocity  profile 
(61:176).  However,  theoretical  and  experimental  studies  of  pulse 
propagation  in  pipe  flow  show  that  local  maxima  or  minima,  other  than 
at  the  centerline,  occur  very  close  to  the  wall  where  the  molecular 
viscosity  is  dominant  (63:83;  70:594-595;  74:867).  As  a  result,  the 
gradient  transport  formulation  should  provide  reasonable  accuracy  for 
the  simple  geometry  and  rigid,  nonporous  walls  of  the  transmission 
lines  in  this  study. 

Turbulence  models  based  on  differential  equations,  such  as  the 
method  of  Bradshaw  et  al.  (6)  and  two-equation  models,  require  the 
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solution  of  differential  equations  in  addition  to  the  governing 
equations  to  determine  a  value  for  the  eddy  viscosity.  This  translates 
into  a  significant  amount  of  computer  time  in  addition  to  that  required 
to  solve  the  governing  equations.  Schetz  (61:208)shows  that 
formulations  based  on  turbulent  kinetic  energy  such  as  that  of  Bradshaw 
et  al.  (6)  are  equivalent  to  the  algebraic  mixing  length  or  eddy 
viscosity  methods  for  steady  flow  in  the  wall  region.  Since  the 
predominant  dynamic  shear  region  is  always  very  close  to  the  wall  in 
transient  flow  in  transmission  lines  (63:84;  70:595;  74:866),  use  of 
the  much  simpler  algebraic  mixing  length/eddy  viscosity  turbulence 
models  for  this  study  is  Justified  on  this  basis.  Anderson,  Tannehill, 
and  Fletcher  (2:226)  point  out  that  algebraic  models  have  accumulated 
an  impressive  record  of  good  performance  for  simple  flows  and  can  be 
modified  for  simple  compressible  flow  applications.  The  steady-flow 
algebraic  eddy  viscosity  model  developed  by  Van  Driest  was  modified  by 
Tellonis  (65)  to  account  for  transient  flow.  Lack  of  data  on  transient 
turbulent  flow  makes  development  of  criteria  for  selecting  a  turbulent 
viscosity  model  impossible.  As  a  result,  the  algebraic  model  described 
by  Tellonis  will  be  used  to  calculate  turbulent  viscosity  in  this 
study. 

The  steady-state  turbulent  viscosity  model  of  Van  Driest  is  used 
as  a  point  of  departure  for  developing  a  model  capable  of  describing 
eddy  viscosity  in  transient  turbulent  flow.  The  steady-state  model  of 
Van  Driest  is  shown  below  (73:477): 


2 
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1  2 

2 

au 

■ 

=  P  < 

1  -  exp 

26  V 

V  / 

y 

ay 
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where  u^= 

This  eddy  viscosity  model  is  valid  for  steady  flow  but  can  be 
modified  to  include  transient  changes  in  the  local  axial  pressure 
gradient  and  velocity.  Cebeci  accounts  for  unsteady  effects  by  using 
T  instead  of  t  in  the  calculation  of  u  (12:2153).  t  is  the  shear 

•  w  *  s 

stress  evaluated  at  the  edge  of  the  viscous  sublayer  rather  than  at  the 
wall.  is  obtained  by  solving  Eq  (17)  which  is  the  momentum  equation 
for  steady  flow  in  the  immediate  vicinity  of  the  wall.  Cebeci 
justifies  this  modification  on  the  physical  grounds  that  the  random 
velocity  fluctuations,  damped  according  to  the  mechanism  of  the  Stokes 
layer,  do  not  extend  all  the  way  to  the  wall  but  disappear  at  the 
viscous  sublayer: 


dr  dp 

S 

dy  dx 


(17) 


Telionis  recommends  the  unsteady  velocity  term  be 
approximation  of  the  momentum  equation  by  modifying  Eq 


included  in  the 
(17)  (66:233; 


77:470): 


du 

dp 

dr 

P  —  = 

—  -  -f  - — 

dt 

dx 

dy 

is  found  by  integrating  from  the  wall  to  to  the  edge  of  the  viscous 
sublayer: 


^  _  1_ 

ax  +(r  -  a)  —  +  p 

tJ  ’  u 

... 

—  r  dr 

S 

r 

W  S  o  ^ 

dx 

dt 

s 

- 

r 

s  J 
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A  more  detailed  development  of  the  eddy  viscosity  model  Is  presented  In 
Appendix  D. 

White  (73:471)  describes  a  three-layer  concept  for  boundary 
layers.  There  Is  an  Inner  layer  where  viscous  shear  dominates,  an 
outer  layer  where  turbulent  shear  dominates  and  an  overlap  layer  where 
both  types  of  shear  are  Important.  The  eddy  viscosity  model  of  Clauser 
Is  usually  assumed  to  hold  In  the  outer  layer  for  flat  plate  flow 
(73:478): 


c  =  p  K  U  5*  (20) 

2  c 

0 

where  =  0.0168  and  5  Is  the  displacement  thickness.  White  (73:473) 
points  out  that  the  outer  layer  Is  commonly  neglected  In  pipe  flow. 

The  outer  layer  has  little  effect  on  the  analysis  and  the  lack  of 
Intermlttence  In  duct  flow  produces  greatly  Increased  turbulence  at  the 
centerline  as  compared  to  flow  over  a  flat  plate.  As  a  result,  Eq  (19) 
will  be  used  to  calculate  eddy  viscosity  from  the  wall  through  the 
Inner  and  overlap  layers  to  the  centerline  of  the  trauismlsslon  line. 

No  special  expression  will  be  used  for  the  outer  layer. 

The  modification  to  the  eddy  viscosity  model  of  Van  Driest 
suggested  by  Telionis  (65)  was  developed  primarily  for  low-frequency 
oscillating  flow.  Pressure  disturbances  examined  with  the 
finite-difference  model  in  this  study  are  not  periodic  and  may  contain 
frequency  components  significantly  higher  than  the  10  -  50  Hz  range  for 
which  the  Telionis  model  was  intended. 

Studies  of  pulsating  flow  make  reference  to  a  frozen  viscosity  in 
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which  the  frequency  of  the  signal  is  much  higher  than  the  mean  turnover 
time  of  the  turbulent  eddies  (11:682,  54:211).  As  a  result,  the  eddies 
do  not  have  time  to  reach  a  new  equilibrium  and  the  change  in  eddy 
viscosity  is  quite  small.  Brown  characterizes  this  frozen  viscosity  by 
a  relaxation  time  constant,  which  is  the  time  it  takes  a  fluid  particle 
near  the  wall  to  travel  a  distance  over  which  its  deviation  of 
turbulence  level  from  its  equilibrium  level  is  reduced  to  e  ^  of  its 
initial  value  (11:686).  This  time  constant  is  somewhat  arbitrarily 
assumed  to  equal  the  time  it  takes  a  particle  to  travel  110  times  the 
boundary  layer  momentum  thickness.  Brown  points  out  that  the  lag  in 
the  establishment  of  turbulence  is  more  important  for  most  problems 
than  the  change  in  the  velocity  profile  due  to  the  dynamic  effects 
(11:686).  He  further  states  that  the  lag  becomes  vastly  more  important 
for  large  disturbances  than  for  small  disturbances. 

The  effect  of  a  lag  in  the  establishment  of  turbulence  in  the 
solution  of  the  equations  developed  in  Chapter  II  could  be  incorporaL’d 
by  treating  the  changes  in  eddy  viscosity  like  the  response  of  a 
first-order  system  where  the  change  in  eddy  viscosity  lags  the  change 
in  velocity  by  a  time  constant  t  Incorporating  a  first-order  lag 

C 

function  into  the  unsteady,  eddy-viscosity  model  is  straightforward. 

Q 

The  ordinary  differential  equation,  (t  =Kp  ,is 

c  ot  lag  T 

solved  numerically  for  with  supplied  by  the  eddy-viscosity 

model.  ^  is  then  used  in  place  of  n  in  Eqs  (5)  through  (7). 

Ia9  ^ 

First-order  behavior  of  eddy  viscosity  could  be  expected  to 
significantly  Influence  the  waveform  of  the  signal.  In  fact,  Brown 
made  an  analogous  modification  to  his  frequency  dependent  model  to 
account  for  eddy  viscosity  in  turbulent  periodic  flow.  The 
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steady-state,  eddy-viscosity  profile  is  multiplied  by  an 
empirically-determined  complex  function  of  frequency  (10). 
Unfortunately,  values  of  x  and  K  could  not  be  calculated  for  use  in 

C 

this  study,  since  no  data  on  unsteady  turbulence  resulting  from  pulses 
or  arbitrary  disturbances  was  found.  As  a  result,  no  provision  for 
relaucation  of  the  eddy  viscosity  will  be  made  in  this  study  although  it 
is  recognized  that  the  effects  of  relaxation  can  be  important. 


Relationships  among  the  various  coefficients  in  the  polynomials 
used  to  approximate  the  velocity  profiles  were  derived  in  Chapter  II. 
These  relationships  provide  a  means  to  deduce  an  initial  steady-state 
velocity  profile  of  am  established  meam  throughflow  without  carrying 
out  the  numerical  solution  of  the  governing  equations  to  steady  state. 

The  mean  flow  is  caused  by  a  difference  in  pressure  between 
opposite  ends  of  the  transmission  line  at  time  t  <  0.  With  the  initial 
value  of  the  pressure  known  at  each  end  of  the  line,  the  steady-state 
average  velocity  for  compressible  flow  in  the  line  is  found  from  the 
one-dimensional  equations  governing  isothermal,  compressible  flow 
through  a  pipe  (4:100): 


The  subscript  0  denotes  the  end  of  the  line  with  the  higher  pressure 
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while  the  subscript  L  denotes  the  low-pressure  end.  The  average 
velocity  of  the  mean  flow  considered  in  this  study  will  always  be 
significantly  less  than  the  speed  of  sound.  Isothermal  behavior  of  the 
fluid  is  assumed  for  the  mean  flow  since  the  mean  velocity  is  low  and 
the  length  of  the  transmission  line  is  much  greater  than  its  radius. 

As  a  result,  the  isothermal  equation  is  used  to  predict  the  mean 
velocity.  The  Darcy-Weisbach  friction  factor,  f,  is  calculated  as  a 
function  of  Reynolds  number  as  follows: 


Laminar  Flow 


(22) 


Turbulent  Flow 


1.14  -  2  log 

21.25 

rO.9 

• 

- 

(23) 


UqD  p 

where  R  =  -  (4:50).  Equation  (21)  cannot  be  solved  for  u  in 

e  fX  0 

closed  form  since  f  is  a  function  of  the  unknown  u  and  an  iterative 

0 

solution  of  Eqs  (21)  and  (22)  or  (21)  and  (23)  for  launinar  or  turbulent 
throughflow,  respectively,  is  required.  Newton’s  method  is  used  to 
solve  these  equation  pairs  for  u^: 


n+l 

U 

0 


du 

0 


The  function  g 


is  obtained  from  Eq  (21): 


(24) 
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2 

U 

0 


(25) 


The  first  derivative  of  Eq  (25)  with  respect  to  is 


du 

0 


-2 


L  _3f 
D  au 


2  u 

0 


(26) 


The  derivative  of  the  friction  factor  is  easily  found  for  both  laminar 
and  turbulent  flow: 


Laminar 


af 

au  “ 

0 


L 

u 

0 


Turbulent 


ar 

du 

0 


3/2 

-  3.6  f  log(e) 


u 

0 


(27) 


(28) 


Iteration  is  used  to  find  u^  from  Eqs  (24)  through  (28). 

Once  u  has  been  determined  from  the  known  values  of  p  and  p  , 
the  pressures  at  each  grid  point  along  the  interior  of  the  line  are 
found.  Equation  (21)  must,  again,  be  solved  by  iteration  but  this  time 
for  pressures  at  each  interior  point  x^,  with  u^  and  f  known.  The 
functions  needed  for  Newton’s  method  in  this  case  are 
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(29) 


(30) 


Values  of  average  velocity  at  each  grid  point  along  the  line  are  found 
from  the  corresponding  pressure  value  at  that  point: 


u  =  u  — 

X  op 

J  X 


J 


(31) 


The  average  meain  velocity  and  pressure  fields  for  a  steady,  mean 
flow  are  used  to  find  steady-state  values  of  the  averaged  velocity 
function  in  each  flow  interval.  The  same  eddy  viscosity  model 
developed  earlier  in  this  section  must  be  used  to  obtain  steady-state 
values  of  the  averaged  velocities  at  the  required  locations  along  tha 
'ine.  The  use  of  a  empirical  expressions  for  the  steady-state  velocity 
profile  produces  an  initial  velocity  field  significantly  different  from 
what  would  be  obtained  from  the  numerical  solution  of  the  governing 
equations  carried  out  until  steady  state  is  reached.  The  use  of  this 
incompatible  field  as  a  set  of  initial  conditions  causes  the  numerical 
solution  to  change  significantly  with  time  in  the  absence  of  an  input 
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signal.  This  drift,  vinder  steady-state  conditions,  would  show  up  as  an 


error  in  the  numerical  solution  when  a  signal  is  present. 

The  relations  for  obtaining  the  coefficients  of  the  polynomials 
described  in  Chapter  II  provide  the  means  to  obtain  an  initial  velocity 
field  compatible  with  Eqs  (5)  and  (6)  or  (10)  amd  (11).  The  key  to 
this  procedure  is  recognizing  that  the  value  of  the  forcing  function  F 

n 

is  the  same  at  the  boundary  of  each  flow  interval  for  steady  flow.  The 
expression  for  F  developed  in  Chapter  II  is  shown  below: 


Values  of  F  that  differ  from  boundary  to  boundary  would  signify  a  net 

n 

change  in  momentum  inside  each  flow  interval  which  is  characteristic  of 
aui  evolving  flow.  Equal  values  of  F  at  each  interval  boundary  also 

n 

imply  that  the  value  of  F  at  each  boundary  is  equal  to  the  value  of  F 

n  M 

at  the  wall.  The  value  for  F  in  the  case  of  an  annular  section  is 

n 

zero  for  steady  flow  when  calculated  as  shown  in  Eq  (32),  since  the 
momentum  entering  the  interval  balances  that  leaving.  Figure  2  shows 
the  flow  cross  section  partitioned  into  both  annular  and  cylindrical 
regions. 

Equation  (32)  may  be  written  at  each  cylindrical  boundary  as 


2  c 


rp  r 


n 


r  1  I I 

r=r  I  dr  J  r  =  r 
n  '  n 


(33) 


Note  that  this  is  Eq  (32)  with  r  equal  to  zero.  For  F  in  Eq  (32) 

n-1  n 

to  be  zero  for  the  annular  sections,  the  value  of  F  in  Eq  (33)  must  be 

n 
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Fig.  2.  Annular  and  Cylindrical  Flow  Partitions 


the  same  at  each  interval  boundary.  The  value  of  F^  at  the  wall  is 
determined  from  the  value  of  the  mean  velocity  (4:42): 


T 

W 


pu^f 

"8“ 


(34) 


This  equation  Is  easily  solved  for 


3r  J  r=a 


to  produce 


3u]  ,  ^ 

dr  j  rsra  8(1 


(35) 


The  result  of  Eq  (35)  is  used  in  Eq  (33)  to  calculate  the  value  of  F^  . 
At  the  wall,  (i^  is  always  equal  to  p.  At  the  interval  boundaries  away 
from  the  wall,  however,  ,  as  described  earlier  in  this  chapter, 
depends  on  the  velocity  gradient  which  must  be  found  by 

or  J  r  =  r 
'  n 

iteration.  A  first  approximation  of  ^  is  obtained  assuming  p^=  p 

■'  n 

aind  solving  Eq  (33): 
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F  r  P  r 
n  n 


r  2  c  (fl  ) 
n  T  r 


r 

n 


(36) 


This  result  is  then  substituted  into  Eq  (16)  to  update  The  new 

value  of  4^  is  substituted  into  Eq  (36)  to  produce  a  better 
approximation  of  ||p|  This  process  is  repeated  until  changes  in 

'  ■'  r=r 

n 

successive  values  of  j^j  are  sufficiently  small  to  indicate 

'■  •'  r  =  r 

n 

convergence.  This  iteration  is  performed  for  all  Interval  boundaries 
excluding  the  boundary  at  the  wall.  The  velocity  gradients  for  a 
laminar  mean  flow  are  obtained  in  a  straightforward  manner  without 
iteration,  since  4  is  known. 

The  values  of  the  averaged  velocity  in  each  flow  interval  can  be 
found  from  the  velocity  gradients  at  each  interval  boundary  with  the 
aid  of  the  recursion  relationships  for  the  coefficients  of  the 
piecewise  defined  polynomials  developed  in  Chapter  II.  The  derivative 
of  Eq  (15)  evaluated  at  the  interval  boundary  provides  an  expression 
for  the  velocity  gradients: 


au 

ar 

J 


r 

n 


/ 

p(x, r  , t )  =  25  c 
1  12 

*  2 
p(x,r  ,t)  =  b+25c+35d 

n  n  n  n  n  n 


n  =  1 

n  =  2. •  ■  •  ,M 


(37) 


where  5  =  r  -  r  Recall  that  b  can  be  expressed  in  terms  of  c 

n  n  n- 1  n  2 

through  c  and  d  as  shown  in  Appendix  C.  Substitution  of  the 

M  M 

recursion  relations  into  Eq  (37)  for  M  =  4  produces  Eq  (38).  Equation 
(38)  is  easily  solved  for  c^,  c^,  c^  and  d^.  The  resulting  values  are 
used  with  Eq  (9)  to  produce  values  of  u  for  steady  flow.  An 
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expression  similar  to  Eq  (38)  can  be  generalized  for  any  number  of  flow 
Intervals: 


[  1^]  1  [  25 

arj  1  1 

|H]  (2S  +5  ) 

drjz  1  2 

1^1  (26  +5  ) 

arJ  3  12 

1^1  (28  +8  ) 

8rj4j  L  1  2 


(38) 
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IV.  Numerical  Scheme 


Two-Step  Lax-Wendroff  Method 

Equations  (5)  through  (7)  and  (10)  through  (12),  with  the 
appropriate  boundary  conditions,  must  be  solved  simultaneously  to 
obtain  the  transient  response  using  a  suitable  numerical  method.  Tsao 
chose  the  two-step,  Lax-Wendroff  numerical  approximation  which  uses 
centered  differences  in  the  spatial  coordinate  to  solve  a  linearized 
form  of  the  damped-wave  equations  (70:590).  Tsao  notes  that  his 
results  for  the  transient  response  of  a  blocked  line  to  a  step  input 
compared  well  with  those  of  Brown  (7;  70:592).  The  Lax-Wendroff  scheme 
works  well  for  the  linearized  damped-wave  equations  and  provides  formal 
second-order  accuracy.  However,  it  does  not  perform  as  well  as  other 
methods  with  respect  to  numerical  diffusion  when  the  local  speed  of 
sound  is  not  constant. 

The  linear,  damped-wave  equations,  Eqs  (10)  and  (12),  written  in 
matrix  form  without  damping  terms  are  shown  below: 


At/2 
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(40) 


n»l 


■  q. 


n*\/2  n*\/2 

q  -  q 
^J+l/2  ^J-1/2 


(Step  2) 


=  0 


(41) 


At 


Ax 


The  accuracy  of  finite-difference  algorithms  can  be  compared  by 
finding  the  modified  equation  and  determining  the  amplification  factor 
for  each  algorithm.  The  modified  equation  is  the  partial  differential 
equation  that  is  actually  solved  when  a  finite-difference  method  is 
applied  to  the  original  partial  differential  equation.  The 
eunplif ication  factor  shows  the  dissipation  and  dispersion  errors 
resulting  from  the  application  of  the  finite-difference  algorithm. 

A  scalar  form  of  the  modified  equation  more  clearly  illustrates 
the  behavior  of  the  numerical  algorithm.  Equation  (39)  is  written  as 

T  q^  +  T  A  T"^T  =  0  (42) 


The  source  of  the  matrix  T  will  be  explained  during  discussion  of 
alternative  finite-difference  methods  later  in  this  chapter.  A  has  two 
linearly  independent  eigenvectors  making  it  similar  to  the  diagonal 
matrix  [.''.^1  containing  the  eigenvalues  of  A  so  that  T  A  T"^=  [A^]. 
Equation  (42)  produces  a  system  of  two  uncoupled,  scalar  equations. 

The  expression 


aq 

at 


+ 


p 

c 


(43) 


is  written  from  Eq  (42)  to  determine  the  most  restrictive  conditions 
regarding  accuracy  and  stability  for  algorithms  applied  to  Eq  (39).  p 

C 


38 


is  the  spectral  radius  of  A  (18:228)  and  is  the  value  of  the  modulus  of 
the  eigenvalue  with  the  largest  modulus.  Application  of  the 
Lax-Wendroff  algorithm  to  Eq  (43)  produces 


n+l  n 
q  =  q  - 


u  iji  c  g  n 
2 


The  modified  equation  is  obtained  from  Eq  (44)  by  substituting  the 
Taylor-series  expansions  for  ^  8®^ 


At 

r 


(At)^  (At)^  ^  ,  ,2  At 

r  24  r 


+  (n  At(Ax)^_ 

2  ^xxx  24  ^xxxx 


(45) 


All  derivatives  with  respect  to  time  are  eliminated  from  Eq  (45)  using 
Eq  (45)  itself,  to  produce  the  scalar  form  of  the  modified  equation 
(2:101): 


(Ax)‘ 


(Ax)' 


(1  -  i;  )  q  -p- 


XXX  c 


1^(1  -  )  q  + 

xxxx 


(46) 


where  v  is  the  Courant  number  (21:212),  which  is  written  as  p  7—  .  If 

Ax 

At  and  Ax  are  chosen  so  that  i'  =  1  ,  the  algorithm  described  by  Eqs 

(41)  amd  (42)  provides  the  exact  solution  to  Eq  (40)  since  the 
right-hand  side  of  Eq  (43)  is  zero.  As  a  result,  the  Lax-Wendroff 
algorithm  satisfies  the  shift  condition  (2:92).  At  v  =  1 ,  this  scheme 
is  equivalent  to  solving  Eq  (43)  exactly,  using  the  method  of 
characteristics. 
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The  extent  of  the  error  resulting  from  the  application  of  the 
Lauc-Wendroff ,  finite-difference  method  is  shown  by  the  amplification 
factor.  For  a  computational  domain  of  length  L  in  the  x-direction,  von 
Neumann  stability  analysis  (27:525)  assumes  the  roundoff  error  is 
written  as 


JaAt  J(i3) 

e(x.t)  =  j;  e  e  (47) 

J 

where  ^  =  k  Ax  ,  k  =  ,  and  j  =  0,l,*»*N  J  =  1 

denotes  the  fundamental  frequency  whose  period  is  2L  and  M  is  the 
number  of  Ax  Intervals  contained  in  length  L.  Equation  (47)  must 
satisfy  the  finite-difference  algorithm  shown  in  Eq  (44).  A  single 
term  of  Eq  (47)  is  written  as  (2:72)  and  is  used  to 

generate  the  amplification  factor  G  for  Eq  (44)  shown  below  (2:102): 

G  =  l-v^(  l-cos0)-tvsin/3  (48) 

Note  that  each  value  of  0  represents  a  Fourier  mode  in  von  Neumann 
stability  analysis. 

Equation  (48)  can  be  written  as  G  =  1G|  e^  where  |G|  is  the 

magnitude  c  f  C  ^  is  the  u.r  For  Eq  (44)  be  stable, 

|G|<1  for  all  values  of  ^  from  0  to  2n.  Values  of  |GI  greater  than  one 

will  allow  roundoff  error  to  grow  without  bound.  Figures  3  euid  4  show 

plots  of  |G|  and  the  relative  phase  error.  The  total  dissipation  from 

the  application  of  the  two-step,  Lax-Wendroff  method  may  be  found  from 

(1  -  |G|'')A  where  A  is  the  initial  amplitude  of  the  wave  and  n  is  the 
0  0 
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number  of  time  steps.  Similarly,  if  represents  the  analytical  phase 
of  the  left-hand  side  of  Eq  (46),  the  total  dispersion  error  of  the 
Lax-Wendroff  method  after  n  time  steps  is  found  from  -  <p) .  The 

total  dispersion  error  0  -  ^  is  the  difference  between  the  phase  of 
the  finite-difference  solution  and  the  phase  of  the  exact  solution  of 
the  differential  equation.  The  relative  phase-shift  error  and  the 

amplification  factor  modulus  |G|  will  be  used  to  characterize 
dispersion  error  and  the  dissipation  error  respectively,  for  a  single 
time  step. 

Plots  of  |G|  and  0/^  are  used  to  compare  the  performance  of  other 
finite-difference  methods  to  the  performance  of  the  Lax-Wendroff 
method.  Figure  3  shows  that  the  Lax-Wendroff  method  exhibits  a  lagging 
phase  error  for  most  values  of  p.  A  relative  phase  error  greater  than 
one  indicates  the  numerical  solution  moves  faster  than,  or  or  leads  the 
exact  solution.  For  a  relative  phase  errors  less  than  one,  the 
opposite  is  true. 

The  use  of  finite-difference  methods  on  long  transmission  lines 
requires  many  time  steps  to  generate  a  solution,  allowing  dissipation 
amd  dispersion  errors  to  accumulate.  Additionally,  when  the  local 
speed  of  sound  varies  with  signal  amplitude,  the  Courant  number  cannot 
be  chosen  to  be  1  everywhere  along  the  line.  The  value  of  v  may  be  1 
at  the  location  of  highest  pressure  but  will  be  less  than  1  at  other 
locations  in  the  transmission  line.  To  maintain  accuracy,  the  total 
dissipation  and  dispersion  errors  must  be  small  for  a  range  of  Courant 
numbers  less  than  or  equal  to  1.  Methods  that  produce  the  smallest 
numerical  errors  when  v  <  1  should  be  chosen  as  alternatives  to  the 
Lax-Wendroff  method  for  the  damped-wave  equations.  Maintenance  of 
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accuracy  in  the  nonlinear  case  is  what  motivates  the  use  of  an 
alternative  algorithm  to  the  Lax-Wendroff  method. 

In  a  liquid,  where  the  local  speed  of  sound  can  be  assumed 
constant,  the  method  of  characteristics  yields  an  almost  exact  solution 
of  the  governing  equations.  However,  if  the  wavespeed  is  strongly 
dependent  on  the  local  pressure  or  temperature,  as  it  is  in  a  gas, 
interpolation  must  be  used  with  the  method  of  characteristics  when 
applied  to  a  fixed  grid.  Traditional  application  of  the  method  of 
characteristics  uses  a  free-running  grid  (1).  The  position  of  the 
point  on  a  fixed  grid  must  be  interpolated  from  points  on  the 
free-running  grid  whose  location  is  not  known  until  after  the 
calculations  are  complete.  Manning  uses  such  a  scheme  on  a  rectangular 
grid  by  approximating  the  characteristic  curves  locally  by  straight 
line  segments  (41).  It  is  this  local  linearization  of  the 
characteristic  curves  and  subsequent  interpolation  of  boundary  points 
that  reduces  the  accuracy  of  this  method  when  applied  to  a  fixed 
rectangular  grid.  This  interpolation  reduces  the  accuracy  of  the 
method  of  characteristics,  making  alternate  numerical  methods  more 
attractive. 

A  number  of  finite-difference  methods  for  solving  the  unsteady, 
Euler  equations  for  compressible  flow  have  been  developed  (15;  20;  45). 
The  finite-difference  method  developed  by  Gabutti  (21)  was  chosen  for 
this  study  due  to  its  small  numerical  dissipation  and  dispersion  errors 
for  the  flow  conditions  under  consideration.  The  basis  for  this 
selection  is  discussed  below. 

The  damped-wave  equations  developed  in  Chapter  II  must  be 
transformed  to  make  use  of  the  finite-difference  methods  developed  for 
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the  unsteady  Euler  equations.  Eqs  (5)  and  (6)  are  written 

q  +  A  q  =  y  (49) 


Ln(p)' 

.  A  = 

■  u  7 

av 

2 

,  and  y  = 

"  H 

a  V 

u 

C 

-  u 

F 

^  n 

1  7 

^  n  J 

where  q  = 

Equation  (49)  can  be  transformed  to  the  required  form,  as 

described  by  Chakravarthy,  by  a  process  similar  to  that  used  to  obtain 

the  compatibility  equations  for  the  method  of  characteristics  (15:2-4). 

Following  Chakravarthy,  the  eigenvalues  of  the  matrix  A  are  found  to  be 

X  =  u  ±  c  .  A  left  eigenvector  is  associated  with  each  eigenvalue  X. 

11 

For  X  “  u  +  c 

n 


■  -c  7  ' 

'  1 

2 

=  0  producing  1  ^  -  a. 

C 

-  -c 

1 

_  r 

c 

(50) 


where  a  is  an  arbitrary  multiplier  which  arises  because  the  length  of 
the  eigenvector  is  not  unique.  Similarly,  for  X  =  u  -  c  ,  the  left 

n 

eigenvector  is  found: 


’  C  7  ' 

'  -1 

2 

=  0  producing  ^2~  ^ 

C 

7 

-  C 

7 

c 

where  0  is  also  an  arbitrary  multiplier.  L  is  a  matrix  whose  columns 
are  the  two  left  eigenvectors  and  with  the  multipliers  a  and  |3 
set  equal  to  one. 
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L  = 


1  -1 

1  1 

c  c 


-1 


and  T"^= 


(52) 


Eq  f49)  is  multiplied  by  T  to  get: 


1  i-* 


-1  I 
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Ln(p) 


1  I 


•1  I 
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Ln(p) 
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1  n 


■1  I 


H 


(53) 


Two  equations  result  from  Eq  (53): 


a(Ln[p] )  7  du 


-  + - “  +  (u  -t-c) - 

at  c  at  "  ax 


a(Ln[p))  y  au  7 

-  +  -(u  +c)— "  =  H  +  -  F 

n  ^  n  n 

c  ax  c 


(54) 


a(Ln(p])  7  au  a(Ln(p])  y  au  7 

-  + - ”  -  (u  -c) -  +  -(u  -c) — ”  =  -H  +  -  F 

at  c  at  "  ax  c  "  ax  "o'* 


(55) 


These  equations  are  used  to  Incorporate  the  boundary  conditions  at  each 
end  of  the  transmission  line  into  the  algorithm  for  the  interior 
points,  u  is  substituted  for  u  in  Eqs  (54)  and  (55)  when  solving 

av  n 

for  Ln(p). 

An  alternate  form  of  the  governing  equations  for  the  internal 
points  is  needed.  For  each  eigenvalue  there  is  a  ieft  eigenvector 
that  must  satisfy  -  A  0)  =  0  .  which  is  equivalent  to 

A  =  i^A^  =  A^i^  Following  Chakravarthy ,  Eq  (53)  is  written  as 
(15:2-4) 
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(56) 


r^q  +  [A  ]  r^q  =T‘‘y 

t  A  X 

Where  [A^]  is  a  diagonal  matrix  of  the  eigenvalues  of  A. 

Multiplication  by  T  produces 

q  +  T  (A  1  T'^q  =  f  (57) 

t  Ax 

Equation  (49)  has  been  recovered  since  [A^]  is  similar  to  A  (18:196). 
The  original  matrix  A  is  now  written  as  T  [A^]  T~^  .  The  matrix  [A^] 
is  split  into  its  positive  and  negative  parts  to  get 
[A  ]  =  [A  ]*+  [A  ]  A  is  now  written  as: 

AAA 

A  =  T  [A^]V^+  T  =  A*  +  A'  (58) 

where 

u  +c  0  I 

n 

0  0 


0  0 

0  u  -c 

n 

Equation  (49)  is  becomes 

q  *  A*  q  +  A'  q  =5=  (59) 

t  X  X 
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The  matrices  A*  and  A  in  Eq  (59)  are  obtained  from  Eq  (58): 


r  U  +C 
n 

5j(u  +c) 

ZC  n 

r  U  -C 
n 

-^(U  -c) 

ZC  n 

2 

2 

^(u  +c) 

2y  n 

u  +c 

n  1 

A"  = 

~{u  -c) 

2y  n 

u  -c 

n 

2  J 

2 

Note  that  A*  +  A  =  A.  The  A*and  A  matrices  for  the  linear 
damped-wave  equations  Eqs  (10)  and  (11)  are  determined  as  shown  above: 


c  7 

C  7 

2  2 

2  2 

2 

A"  = 

2 

c  c 

C  C 

2y  2 

- 1 

l(M 

1 

1 

_ 1 

The  linear  matrices  differ  only  by  the  lack  of  u  in  each  matrix 

n 

element. 

Substituting  these  expressions  into  Eq  (59)  and  performing  the 
required  multiplication  produces  a  set  of  simultaneous  partial 
differential  equations: 

a(Ln[p])  u  +c  a(Ln(p])  r  du  u  -c  a(Ln[p] ) 

-  + - +  - (u  +C) -  + - - 

at  2  ax  2c  ax  2  ax 

y  au 

—  (u  -c)— *''=  H  (60) 

T  av  av 

2c  ax 


au  c  a(Ln(p]  )  u  +c  au  c  a(Ln[p]  )  u  -c  au 

—  +  —  (u  +c) -  +  -  —  -  —  (u  -c)  -  -»•  -  —  =  F 

at  2y  "  ax  2  ax  2y  "  ax  2  ax  " 

(61) 

Appropriate  one  sided  differences  must  be  used  to  approximate  the 
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spatial  derivatives.  Derivatives  with  coefficients  containing  are 

approximated  using  backward  differencing,  since  these  coefficients 
result  from  the  A*  matrix  and  represent  right  running  waves. 

Similarly,  derivatives  with  coefficients  containing  u  -c  are 

n 

approximated  using  forward  differencing.  This  is  known  as  upwind 
differencing  and  is  always  in  the  opposite  direction  of  propagation. 

Three  finite-difference  methods  applicable  to  Eqs  (60)  and  (61) 
were  considered.  The  method  with  the  least  numerical  dispersion  and 
dissipation  is  sought.  A  first-order  method  by  Chakravarthy  et  al. , 
known  as  the  split -coefficient,  matrix  method  (SCM),  relies  on  a 
two-point  differencing  scheme  (15).  The  finite-difference 
representation  of  Eqs  (60)  and  (61)  is  written  as 

n-*  1  n  At,.*  _  n  ^  .  n. 

q  =  q  -  T— (A  7  q  +  A  A  q  )  (62) 

where  7  q"  =  q”-  q"  and  A  q"  =  q"  -  q"  This  method  is 
J  J  J-i  J 

first-order  accurate  in  space  and  time,  requiring  relatively  few 
operations  to  calculate  successive  points.  However,  a  very  fine  grid 
is  required  to  achieve  good  accuracy. 

Second-order-accurate  methods  require  fewer  computational  nodes  to 
achieve  a  given  accuracy  than  do  first-order  methods  but  more 
arithmetic  operations  are  required  to  calculate  successive  points.  A 
number  of  second-order  methods  suitable  for  application  to  the 
damped-wave  equations  exist.  One  method  proposed  by  Moretti  uses 
three-point  differences  in  a  predictor-corrector  arrangement  (46).  At 
the  predictor  level 
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At  the  corrector  level 


n*l  1/  n  n'*'!  .  r  \ 

‘*J  “  2^‘*j  ’j 

where  q"  denotes  q"  =  A*q~+  A~q*  The  differencing  is  reversed  from 
that  of  the  predictor  level: 


n-*'l  n+l 

q  -  q 

^J-1 


-2q7^ 


3q 


ni'*’! 

J*1 


-  q 


„n+l 

J+i: 


q„  = 


Ax 


and  q  = 


Ax 


(65) 


This  scheme  is  second-order  accurate,  however  it  does  not  satisfy  the 
shift  condition. 

A  family  of  second-order-difference  schemes  that  do  satisfy  the 
shift  condition  has  been  developed  (21).  The  Gabutti  schemes  use  a 
predictor-corrector  approach,  as  does  the  Moretti  scheme,  but  the 
differencing  scheme  is  changed  and  the  predictor  is  subdivided  into  two 
parts: 


q"*^=  q"-  At(A*q  +  A  q* )  (66) 

J  J  XX 
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(67) 


q  =  A  q  +  A  q 


The  spatial  differences  for  the  first  part  of  the  predictor  are 


q 


X 


Ax 


Ax 


(68) 


For  the  second  part  of  the  predictor  the  spatial  differences  are 


♦ 

q  = 


~  n  _  n  n 
-2q  +  3q  -  q 


_  n  _  n  ^  n 

2q  -  3q  +  q 

J-l  J-2 


q  = 


Ax 


Ax 


(69) 


The  corrector  is  written 


n+l  n  At ,  n  n+1 , 
'j  '  ’j  *  -‘"t  *  "t  > 


(70) 


where  q"^^  is  defined  as 


"♦I  ft*  -  „-  * 

q  =  A  q  +  A  q 

t  XX 


(71) 


and 


n+l  n+l 

q  -  q 
J  J-i 


Ax 


n+l  n+l 

^J  +  1 


Ax 


(72; 
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This  same  scheme  was  earlier  developed  by  Warming  and  Beam  (2:294). 
Equations  (66)  through  (72)  are  combined  to  produce 


n.l  n  ^  n  ^  ^  2  n 

J  J  2  ^  ^  2AX  J 


2Ax  ^ 


(73) 


Scalar  expressions  for  numerical  dispersion,  numerical  dissipation 
and  the  modified  equations  for  the  three  algorithms  described  above, 
obtained  as  discussed  for  the  Lax-Wendroff  method,  are  shown  below  for 
the  wave  equation  with  v  >  0: 


Upstream  Method  (SCM)  (2:90-9?) 

<1,  P  =  P  —(1  -  v)q  -  p  —  — -(2t>^-  3v  +  l)q  + 
t  ex  c  2  XX  c  ^  XXX 


G  =  l-i^(i-cosB)+ii'sin/3 


(74) 

(75j 


Moretti  \  Differencing  (21:212) 

P  =  p-^^^(2  +  i/^)q  -  p  (6+3i/-2i^^)q  +  •  • 

t  c  X  CO  XXX  c  24  xxxx 


G  =  l-i’(l-cosB)^-i^^(l+cosB-4cos^/3+2cos^0)- 


ivsinB[2-cosB-2v(  l-cos(3)  ] 


(76) 


(77) 


Upwind  Method 


(Gabutti,  Warming-Beam) 

( Ax ) f  . 

=  p  —7 — (l-v)(2-i')q  -p 

c  O  XXX  c 


(2: 103-104 

(Ax)^,  ,2 

-s —  (1-v) 


;  21:219) 
(Z-vjq  + 

XXXX 


(78) 


52 


l-2i'j^y+2(l-i')sin‘‘  -J  sin^  -  -iysin/3j^  l+2( l-i')sin”  -  J 


Figures  5  thro«igh  10  show  the  dissipation  and  relative  phase 
errors  for  each  of  the  above  finite-difference  methods.  These  figures 
show  the  error  behavior  for  one  time  step.  Up  to  200  time  steps  are 
required  to  obtain  the  transient  response  of  the  systems  considered  In 
this  study.  The  result  Is  continued  accumulation  of  dissipation  and 
phase  error  at  each  time  step.  Based  on  aunpllf Icatlon  factors 
developed  from  the  scalar  equations  shown  above,  the  algorithm  that 
exhibits  the  least  dissipation  error  and  relative  phase  error  was 
selected  for  use  In  this  study.  An  examination  of  Figures  3  through  10 
Indicates  that  the  Gabuttl  method  produces  the  least  dissipation  and 
dispersion  error.  Further,  Gabuttl  states  that  in  the  case  of  subsonic 
flow,  the  phase  errors  produced  by  the  forward  and  backward  differences 
used  in  this  algorithm  tend  to  balance  out  (21:214). 

Equation  (80)  is  the  complete  amplification  ^actor  for  the  Gabuttl 
method  as  applied  to  Eq  (39)  (21:220): 


G  =  0-  -p-  [(A*-A  )  (l-cos/3)^-i  [A  )  sin8(2-cos8 )  ]  + 

ilX  A 

[ax]  )^]cos8(cos8-l )+i  ^  [(A*)^-(A  )^]sin3(l-cos8)  (80) 


Complete  amplification  factors  for  the  other  split  coefficient  methods 
discussed  above  are  not  shown  but  the  expressions  are  of  similar  form 
and  contain  IK*  and  A”. 


Additionally,  Chaudhry  and  Hussain!  compared  the  performance  of 
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Fig.  7.  Magnitude  of  Amplification  Factor  for  Moretti  \  Differencing 
Method 
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the  Moretti  \-scheme,  Gabutti  scheme,  the  method  of  characteristics, 
and  the  MacCormack  method  for  a  waterhammer  wave  produced  by  an 
instantaneous  valve  closure  in  a  pipeline  (16)-  The  system  of 
equations  used  to  compare  the  methods  is  shown  below: 


H 

t 


^  iA  Q.  =  ° 


(81) 


+  gA  +  RQ|Q|  =  0 


(82) 


The  MacCormack  method  is  a  variation  of  the  two-step,  Lax-Wendroff 
method  and  exhibits  the  same  dissipation  auid  dispersion  behavior  as  the 
Lax-Wendroff  method  (2:103).  A  comparison  of  the  error  between  the 
computed  solution  and  the  exact  solution  was  made  for  the  different 
numerical  schemes  using  Courant  numbers  less  than  one.  These  results 
are  important  to  this  study  since  the  Courant  numbers  will  be  less  than 
one  at  points  in  the  transmission  line  removed  from  points  of  peedc 
pressures  when  the  numerical  method  is  applied  to  a  nonlinear  problem. 
The  Moretti  method  produced  unacceptable  post  shock  oscillations.  In 
fact,  Gabutti  shows  that  the  Moretti  method  is  unstable  for  v  >  2/3  , 
so  this  result  is  not  surprising  (21:213).  TTie  performance  of  the 
Gabutti  method,  the  MacCormack  method,  and  the  method  of 
characteristics  compared  well  with  the  exact  solution  at  Courant 
numbers  less  than  one.  This  compai ison  of  methods  did  not  address 
errors  i 'suiting  from  numerical  dispersion. 

The  Chaudhry  and  Hussain!  comparison  of  numerical  methods,  in 
conjunction  with  the  numerical  performance  of  the  three  methods, 
discussed  earlier,  supports  the  Gabutti  method  as  the  bes^"  alternative 
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V.  Lumped  Parameter  Representation  of  Transmission  Line 


Lvunped-Parameter  Models  of  Fluid  Transmission  Lines 


The  components  of  a  finite-difference  model  for  fluid  transmission 
lines  were  developed  in  the  previous  chapters.  While  this  model  Is 
very  useful  for  understanding  transient  behavior  of  single  transmission 
lines,  application  to  fluid  systems  containing  many  components  and 
other  transmission  lines  requires  computer  capability  normally  not 
available  to  the  design  engineer.  Practical  fluid  system  design 
requires  an  engineering  model  that  accounts  for  transmission  line 
transient  behavior  in  the  total  system  response  while  requiring  only 
modest  computer  resources. 

A  transmission  line  is  usually  part  of  a  complex  system  consisting 
of  devices  represented  in  terras  of  time  only.  Unlike  the 
distributed-pars.aieter  model  developed  in  previous  chapters,  the 
lumped-parameier  model  describes  the  dynamic  behavior  of  a  fluid 
transmission  line  in  terms  of  time  only.  A  number  of  lumped-parameter 
approximations  for  the  linear  transmission-line  model  have  been 
developed.  One  particularly  useful  lumped-parameter  model  is  based  on 
modal  superposition  (31;  40).  This  formulation  results  in  a  system  of 
equations  compatible  with  the  state-space  method  commonly  used  in 
calculating  ♦he  transient  response  of  control  systems. 

Modal  superposition,  which  is  based  on  classical  mechanics, 
approximates  the  total  transmission-line  response  by  superposing  the 
transmission-line  response  for  each  of  its  normal  modes.  A  system  with 
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n  degrees  of  freedom  will  have  n  natural  frequencies.  Under  free 
vibration  at  one  of  its  natural  frequencies,  a  definite  relationship 
exists  among  the  amplitudes  of  each  of  the  n  coordinates  corresponding 
to  each  of  the  n  frequencies  for  a  mechanical  system  such  as  a 
prismatic  beam.  This  relationship  is  known  as  a  normal  mode  of  the 
beam  (62:16-17). 

Modal  superposition  describes  the  dynamic  properties  of  the 
transmission  line  in  terms  of  its  modes  of  "vibration".  The  modes  of 
vibration  can  only  be  defined  for  linear,  dynamic,  systems  but 
nonlinear  effects  may  be  empirically  incorporated  into  modal 
superposition  by  adjustment  of  the  model  coefficients  (47).  In  the 
ainalysis  of  structural  transients,  the  description  of  a  structure’s 
dynamic  properties  by  a  normal-mode  model  is  computationally  efficient 
since  it  Involves  relatively  few  parameters.  The  dynamic  properties  of 
a  fluid  transmission  line  may  also  be  defined  by  specifying  its  modes 
of  "vibration."  Vfhen  the  dynamic  properties  of  a  structure  or  a 
transmission  line  have  been  determined,  they  can  be  used  to  predict 
transient  responses  to  an  arbitrary  excitation. 

Simplified  Equations  for  Solution  by  Normal  Modes 


The  same  principles  used  to  model  the  longitudinal  vibration  of 
prismatic  beams  in  classical  dynamics  may  be  applied  to  a  fluid 
transmission  line.  Tsao  (70)  developed  a  set  of  linearized  damped-wave 
equations  which  provide  a  connection  between  the  equations  developed  in 
Chapter  II  and  mod"’^  superposition: 
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+  ^  =  2ga(y-l)  ^  £1] 

rp  3t  ax  ~  yp  T  arjr=a 
*^0  •' 


^  ^  A  ap  _  2na  aul 

Jt  P  ax  ~  P  arjr=a 

^0  0  ■' 


(84) 


Equation  (83)  as  originally  shown  by  Tsao  is  incorrect.  The  right-hand 

2Ka(y-l)T  kRT 

side  of  Tsao’ s  Eq  (19)  becomes  -  when  -  is  rewritten  as  -  . 

jp  <r  c  p 

A  check  of  Tsao’s  development  of  the  equation  shows  the  factor  i  was 

omitted  from  the  second  term  of  Tsao’s  Eq  (9).  Equations  (S3)  and  (84) 

are  in  terms  of  pressure  and  volumetric  flow.  The  dissipation  terms  on 

the  right  hand  side  of  each  equation  present  a  problem  to  modal 

superposition  since  the  information  necessary  to  evaluate  the  velocity 

and  temperature  gradients  is  not  available.  The  left-hand  sides  of  the 

equations  pose  no  problem. 

The  perfect  gas  law  p^  =  expression  for  the  speed  of 

2 

sound  c  =  yRT,  where  R  is  the  gas  constamt,  are  used  to  rearrange 
Eqs(83)  and  (84).  Neglecting  the  heat  transfer  term  and  assuming 
Hagen-Poiseuille  flow,  Eqs  (83)  and  (84)  can  be  written  as  shown  below 
(40:459;  61.20): 


+  i  ^  = 

A  ax  ^  p  at 


(85) 


p  ax  A  at  A 


0 


(86) 


Following  Lebrun,  the  dissipation  resulting  from  friction  is  expressed 
as  FQ  where  F  is  a  loss  coefficient  calculated  as  F  =  8u/a^  for 
Hagen-Poiseuille  flow  (40:459).  Although  heat  transfer  is  important 
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for  transient,  compressible  flow,  no  simple  expression  for  heat 
transfer  in  terms  of  Q  alone  can  be  developed  for  the  right-hand  side 
of  Eq  (83).  A  disturbance  must  be  applied  to  each  end  of  the 
transmission  line  to  introduce  pressure  or  flow  signals.  For  the  case 
of  pressure  signals,  a  forcing  function  ♦  for  use  with  Eq  (86)  is 
written  as 


♦  =  i  (p^5(0)  +  P25(L)]  (87) 

where  5  stands  for  the  Dirac  delta  function. 

As  shown  by  Karnopp  (31),  Eqs  (85)  through  (87)  are  combined  to 
form  a  single  second  order  partial  differential  equation  in  terms  of  Q: 

p  a^Q  p  F  dQ  pc  a\  1 

- -  - - -  -  -  (p  5(0)  +  p  5(L)]  (88) 

A  at^  A  at  A  ax  L  ^ 

The  dots  placed  above  variables  denote  differentiation  with  respect  to 
time. 

Modal  Superposition 

Following  the  development  of  Lebrun  (40),  rename  Q  and  p  as  f  and 
p  respectively.  The  symbol  f  denotes  a  flow  variable  while  e  denotes 
an  effort  variable.  These  variables  arc  the  power  variables  associated 
with  bond  graph  modeling  (42:96).  Bond  graphs  are  a  concise,  pictorial 
representation  of  interactive  system  dynamics  from  which  governing 
system  state—pace  equations  can  be  derived  directly  (42:94).  Bond 
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graphs  also  allow  dependent  and  Independent  variables  to  be  determined 
in  an  organized  mauiner.  Lebrun’s  development  of  modal  superposition 
and  proper  application  of  the  model  to  complex  fluid  systems  rely  on 
bond  graphs.  A  detailed  explanation  of  bond  graph  modeling  techniques 
is  provided  by  Kamopp  and  Rosenberg  (32).  Lebrun’s  modal 
representation  of  Eq  (88)  is  obtained  by  expressing  f  as  the  product  of 
mode  shapes,  G^(x),  times  modal  general  coordinates,  (31:439): 

00 

f(x,t)  =  I  G^(x)  Cj(t)  (89) 

1=0 


is  an  independent  coordinate  analogous  to  the  general  coordinate 
used  to  analyze  the  motion  of  a  mechanical  system  by  modal 
superposition.  Substitution  of  Eq'  (89)  into  Eq  (88)  produces 

P  pF  3^  pc^  d^G 

- i  + - L - 1  =  0  (90) 

AC  at  AC  at  AG  5x 


The  forcing  function  on  the  right-hand  side  has  been  temporarily  set  to 
zero,  since  the  normal  modes  are  determined  under  free  vibration 
(32:82).  From  Eq  (85),  =  0  ,  and  =0  ,  in  the  absence  of 

•'  x=0  '  x=L 

the  forcing  function.  These  are  the  boundary  conditions  for  Eq  (90). 
From  Eq(90),  the  separation  of  variables  technique,  detailed  in 
Appendix  E,  produces  a  system  of  equations  with  the  dependence  upon  the 
spatial  variable,  x,  removed: 


pL  3^C  PLF  aC  pc^i^TT^ 

- i  + - - ^  =  e  +  e  (-1)*  (91) 

2A  at  2A  at  2AJ.  '  ^  ^ 
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The  variables  e^  and  e^  denote  pressure  Inputs  at  x  =  0  and  x  =  L, 
respectively.  Integration  over  time  produces  the  following: 

pL  pLF 

- -  *  -  C  =  e  +  e  (92) 

A  at  A  °  ^  ^ 

pL  ac,  pLF  pc^l^ir^t 

- -  +  -  C  +  -  c  dt  =  e  +  e  (-1)  (93) 

2A  at  2A  ‘  2AL  -^0  ‘  ^ 

Karnopp  points  out  that  Eq  (92)  represents  the  line  as  a  single, 
lumped.  Inertia  while  Eq  (93)  represents  a  collection  of  modal 
oscillators  (31:439).  The  following  electrical  analogies  are  drawn  for 
this  system  of  equations  (32:40): 


where  ♦  ,  ♦  ,  and  ♦  symbolize  resistance,  capacitance,  eind  inertance 

R  C  I 

elements  respectively.  Equations  (94)  through  (96)  allow  the  inertia, 
capacitance,  and  resistance  parameters  in  Eqs  (92)  and  (93)  to  be 
identified  by  inspection.  These  parameters  are  listed  in  Table  3  where 
R,  C,  and  I  denote  values  of  $  ,  ♦  and  ♦  ,  respectively.  The  stiffness 

R  C  I 

K  is  defined  as  1/C. 

The  significant  frequency  components  of  the  input  signals  are 

likely  to  be  contained  inside  some  bandwidth.  Let  0  to  w  be  the 

0 

frequency  range  containing  the  significant  frequency  components  of  the 
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(94) 

(95) 

(96) 


Table  3 


Inertia,  Capacitance,  and  Resistance  Parameters  Pressure  Inputs 


I 

C 

R 

pL 

pLF 

l»0 

— 

00 

— 

A 

A 

1=1,2, 

pL 

2AL 

pLF 

2A 

2  2,2 
s  pc  i 

2A 

pressure  signal  imposed  at  the  ends  of  the  transmission  line.  There 
are  an  infinite  number  of  modes  as  shown  in  Table  3,  however,  all  modes 
above  may  be  treated  as  stiffness  controlled  (32:84).  This 
approximation  allows  all  frequency  modes  above  to  be  represented  by 
a  residual  compliance.  The  effects  of  all  modes  above  the  frequency  of 
the  highest  mode  retained  are  summed  up  using  this  residual  compliance. 
An  expression  for  residual  compliance  is  derived  in  Appendix  E: 


C 

oq 


2AL 


2  2 
pc  TT 


J  =  1 
00 


(n^  J) 


=0 


(_i 

E  ^.2 

J  =  1 


( n+  J  ) 


^  TiiTJ) 
J  =  1 
00 

(n+j: 


1  =  1 


(97) 


Equations  (95)  and  (96)  may  be  written  in  general  form  by  defining 
a  new  set  of  generalized  coordinates  (32:84;  40:461): 


q 


1 


p 


1 


(98) 

(99) 
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In  terms  of  these  general  coordinates,  Eqs  (93)  amd  (94)  are  written  as 


IM^2q  +  R|^q  +  Kq  =  (100) 

Symbols  used  in  Eq  (100)  are  defined  in  Table  4.  A  matrix  form  of  Eq 
(96)  is  written  as 


H  It  ?  =  H  q  =  . 


(101) 


where  ^  is  a  vector  with  components  .  Integrating  once  over  time 
produces 


M 


e  dt 


P 


(102) 


Equations  (101)  and  (102)  allow  Eq  (100)  to  be  written  in  a  state-space 
format: 


p  =  IH%  -  K  q  -  R  M’^p  (103) 

e 

q  =  M“^p  (104) 

where  p  and  q  denote  differentiation  of  p  and  q,  respectively,  with 
respect  to  time.  At  each  end  of  the  transmission  line  there  is  an 
applied  pressure  and  a  resultant  flow.  An  output  expression  for  the 
flows  caused  by  the  application  of  the  pressure  signals  is  needed. 

The  bond  graph  in  Figure  11  represents  Eqs  (103)  and  (104)  with  a 
pressure  input  and  a  flow  output  at  each  end  (40:460).  The  variables 
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Table  4 


Definition  of  terms  in  modal  state  equations  (39:461) 
<1  *  "  ■  <I  j  generalized  displacement  vector 

p  =  Pj  * ' '  P  j  generalized  momentum  vector 
e  =  ^e  ^  e  input  effort  vector 

f  =  ^  f  output  flow  vector 


M  = 


K  = 


IR  = 


IH  = 


0 

0 

K 

1 

0 

R 

0 

R 

i 

0 

1  1  1 

1  -1  1 


inertauice  matrix 


stiffness  matrix 

2,2  2 

pc  L  n 


where 


2AL 


R 


resistamce  matrix 


(-1) 


coefficients  of  input  vector 


chosen  as  input  variables  are  clearly  shown  on  the  diagram  and 

determine  the  output  expression.  The  subscript  e  denotes  inputs  whil 

the  subscript  s  denotes  outputs.  The  output  relations  for  f  and  f 

sl  s2 

are  obtained  directly  from  Figure  11  by  performing  the  summations 
indicated  by  the  0-node  at  each  end  o^ 


rv  w 

%  ^ 

TF  TF 


-I  c  (. 

'  ea  I 


0 


— 

s2 


Fig.  11.  Bond  Graph  for  Eqs  (106)  and  (107)  with  Pressure  Inputs  at 
Each  End 


the  transmission  line.  The  result  is  then  integrated  with  respect  to 
time,  dt  and  dt  are  treated  as  displacements 

al  s2 


'  -t  09 

J 

Z  Cj(t)  dt 

1=0 

r*- 

e 

e  1 

t  00 

r  1  (-i)‘?j(t)  dt 
•  ''  1=0 

C  f 

eq  J 

f  dt 

L  J  S2  J 

e 

e2 

Displacement  Sun  of  Displacements  Displacement  Due  to 

at  Ends  of  Line  from  Individual  Modes  Application  of  Force 


(lO: 


Solving  Eq  (105)  for  f  produces  the  output  expression  for  the 

S 

transmission  line  with  a  pressure  applied  to  both  ends 


r 


f  dt  =  W 

S 


C  e 

eq  e 


(lOi 


The  IH  matrix  in  Eq  (106)  results  from  the  factorization  of  the  first 


term  on  the  right-hand  side  of  Eq  (105) 


1  1  •  • 


1 


-t  « 

J 

laO 


(t)  dt 


*  I 

I  Z  (-1)  €.(t)  dt 

**  laO 


1  -1 


An  alternate  output  expression  is  required  for  a  transmission  line 
having  a  pressure  input  at  one  end  and  a  flow  input  at  the  other.  For 
this  case,  the  vectors  e  and  f  now  contain  the  elements  e  ,  e  and 

el  a2 

f^^,  respectively.  Each  vector  has  one  input  and  one  output 
variable.  The  blocked  line  is  an  example  of  this  case.  The  flow  input 
at  the  blocked  end  is  constant  at  zero  while  a  pressure  signal  is 
applied  at  the  open  end.  Figure  12  shows  the  bond  graph  for  this  case. 


Fig.  12.  Bond  Graph  for  Mixed  Pressure  and  Flow  Inputs 


Following  Lebrun,  Eq  (106)  is  first  solved  for  f  to  produce: 

si 


22 


C 


21 


C 


22 


f 

e2 


(108) 
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where  L  H  and  L  H  denote  the  upper  row  and  lower  row  of  matrix  H, 

1  2 

res|)ectlvely.  Equation  (106)  is  then  solved  for  e^^  with  Eq  (108) 
substituted  for  f  in  the  result: 


e 

■  2 


1 


22 


dt  -  L  H 
2 


(109) 


Equations  (103)  and  (104)  with  Eqs  (108)  and  (109)  as  output  relations 
form  the  system  of  simultaneous  equations  governing  the  dynaunics  of  the 
blocked  line  shown  in  Figure  12. 

The  expressions  for  1,  C,  and  R  shown  in  Table  3  are  of  the  same 
form  as  ..hose  contained  in  the  average  friction  model  for  a  circular 
line  with  no  heat  transfer  f36;78): 


I  * 

a 


P 

A 


A 


C  = 

a 


Snfi 


(110) 


(111) 


(112) 


Equations  (91)  and  (93)  are  rewritten  in  terms  of  Eqs  (110)  through 
(112)  and  the  generalized  coordinater  defined  in  Eqs  (98)  and  (99): 


d^q  dq 

L  I  - L  R  — °  =  e  +  e 

*  dt^  dt  ^  ^ 


(113) 
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L  a  q,  L  dq  I  n  1 

-  I  -  +  -  R  — ‘+ - q  =  e  +  e  (-1)‘  (114) 

2  *  at^  2  at  2  L  C  ‘  ^  ^ 

a 

Equations  (113)  and  (114)  describe  the  linear  behavior  of  a  fluid 
transmission  line  assuming  the  parameters  (inertance),  R 
(resistemce) ,  and  (conductance)  are  known.  Similar  equations  may  be 
obtained  for  an  electrical  transmission  line. 

Transmission  line  theory  requires  four  distributed  parameters  to 
characterize  a  line  (37).  Figure  13  shows  am  infinitesimal  section  of 
length  Ax  of  a  uniform  transmission  line. 


lAx  RAx  i(x+Ax,t) 

* 

w  w - ^  NX  - 1 

v(x,t)  < 

< 

> 

—  CAx 

v(x+Ax, t ) 

i(x,t)  l(x+Ax,t) 

Fig.  13.  Equivalent  Circuit  of  an  Infinitesimal  Portion  of  Uniform 
Transmission  Line  in  Time  Domain 


Electric  circuit  theory  produces  two  partial  differential  equations 
from  Figure  13: 


dv ( X , t ) 
dx 


R  i(x,t)  -  I 


5i (x, t ) 

at 


(115: 
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(116) 


dl(x,t)  _  _  p  irfv  t  1  -  C 

where  v  denotes  voltage  and  1  denotes  current.  These  two  expressions 
are  combined  to  form  one  equation  In  the  same  manner  that  Eq  (88)  was 
obtained  from  Eqs  (85)  and  (86): 

a^Kx.t)  fi  G  ai(x,t)  1  a^Kx.t)  r  g 
I  -  +  - R - -  -  Kx.t)  =  0 

at^  c  at  c  ax^  c 

V  I 

(117) 

Equation  (117)  Is  similar  In  form  to  Eq  (88)  with  the  forcing  function 

set  to  0.  In  Eq  (117),  1  replaces  Q,  I  replaces  ^  ^  replaces 

1  oc 

and  ^  replaces  ^  .  The  modal  representation  of  Eq  (117)  with  the 
forcing  function  set  to  zero  Is  obtained  by  the  same  method  used  to 
obtain  Eqs  (92)  and  (93): 

a^q  aq 

L  I  - -  +  L  R  —  =  0  (118) 

at^  at 

L  a^q  L  aq  L  f  1  R  G  ' 

-  I  - -  +  -  R  _Ji  +  - - +  -  q  =  0  (119) 

2  at^  2  at  2  [  C  C  J  “ 

00 

where  q  =  dt  l(x,t)  =  T  G  (x)^  (t)  and  G  (x)  Is  the  mode 

k  k  ,  k  k  ,  k 

k  =  0 

shape.  Comparison  of  these  equations  with  Eqs  (113)  and  (114)  suggests 
a  way  to  include  the  parameter  G  into  the  modal  representation  of  the 
fluid  transmission  line  provided  there  is  some  means  to  calculate  G. 
Note  that  Eqs  (118)  and  (119)  reduce  to  Eqs  (113)  and  (114)  for  G  =  0. 
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An  approximation  of  the  modal  frequencies  for  the  transmission 
line  is  obtained  using  the  expressions  for  I  and  C  shown  in  Table  3 
with  Eqs  (103),  (104),  (108),  and  (109): 


[l  —  L  H 

I  2  2 


0 


(120) 


Detailed  development  of  this  expression  is  shown  in  Appendix  E. 

Solution  of  Eq  (120)  produces  an  approximation  of  the  line’s  modal 
frequencies  for  whatever  number  of  modes  considered.  Note  that  the 
effect  of  damping  has  been  neglected  and  G  has  been  assumed  to  equal 
zero.  More  accurate  expressions  for  I  and  C  produce  better 
approximations  of  the  modal  frequencies  as  well  as  improve  solution 
accuracy  for  the  transient  response  of  the  transmission  line  to  an 
applied  signal. 

The  expressions  for  I,  R,  and  C,  shown  in  Table  3,  do  not  vary 
with  frequency.  It  is  well  known  that  values  of  these  parameters  as 
well  as  values  for  the  conductzuice,  G  ,  are  strong  functions  of 
frequency.  Frequency  dependent  values  for  these  parameters  can  be 
obtained  from  the  frequency-domain  techniques  developed  to  calculate 
accurate  values  for  I,  R,  C  and  G. 

Since  I,  R,  C,  and  G  are  functions  of  frequency,  values  must  be 
calculated  at  the  correct  frequencies  for  the  system  under 
consideration.  One  method  to  obtain  the  correct  modal  frequencies  is 
to  use  Eq  (120)  in  conjunction  with  frequency  dependent  expressions  for 
I,  R,  C,  and  G.  This  method  requires  iteration  and  repeated 
determination  of  the  eigenvalues  from  Eq  (120)  and  ignores  the  effects 
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of  viscous  damping  and  heat  transfer.  A  more  efficient  method 
calculates  the  modal  frequencies  from  the  transmission  line  transfer 
matrix  (36:68): 


P  slnh  n.  +  fZ  /Z  Icosh  FL 

1  I  L  cj 


(121) 


where  Z  is  the  load  impedance,  Z  is  the  characteristic  impedance,  amd 

L  c 

r  is  the  propagation  operator.  The  ratio  calculated  over  a 

range  of  frequencies  containing  the  significant  frequency  components  of 

d  P  /P 
2  1 

the  input  signal.  The  values  of  w  where  both  -  *  0  and 

d(i> 

d^  /P 

2  I 

-  <  0  are  the  modal  frequencies. 

dw^ 

Values  of  Z  ,  Z  ,  and  T  in  Eq  (121)  are  obtained  as  shown  in 
c 

Appendix  A.  Expressions  requiring  evaluation  of  complex  Bessel 
functions  have  been  developed  for  calculating  I,  R,  C,  auid  G  (47:7), 
however,  a  simple,  but  accurate,  approximation  of  these  expressions 
exists.  These  approximate  expressions,  modified  for  use  with  Eqs  (118) 
and  (119),  are  shown  below  (37): 


R(u)  =  -  [DR]  (122) 

.2 


p  8x4 

Ku)  =  -  +  -  [DL]  (123) 

A  A^ 


A  (y  -  1)  P  [DC] 

C(w)  =  -  -t-  - 

^  ^  r  p([DC]^+  [DG]^ 


(124) 
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GluJ 


(y  -  1)  A  [DG] 


y  p| 


[DC]^+ 


[DG] 


1 


(125) 


where 


DR 


3 

8 


*  ih  * 

4  V 


3  i 
8  h 

V 


DI 


1  .  ^  1 
4  %  ■  64  h 


- 1  1  .  ^11 

4  2  T  4  h 

T 


DG 


In 


1 1 
4  h 

T 


(126) 


(127) 

(128) 


(129) 


Expressions  for  and  are  shown  in  Appendix  A. 

Equations  (122)  through  (129)  are  used  to  calculate  values  for  I, 
R,  C,  and  G  at  the  various  values  of  u  as  Eq  (121)  is  evaluated  over 
the  frequency  range  under  consideration.  The  modal  frequencies, 
determined  in  this  manner,  are  then  used  in  Eqs  (122)  through  (129)  to 
calculate  values  of  I,  R,  C,  auid  G  for  use  as  modal  parameters.  The 
values  of  these  four  parameterb  determine  the  tramsient  response  of  the 
transmission  line  to  an  applied  input  signal  when  substituted  into  Eos 
(118)  and  (119).  .  cc r olementary  set  of  equations  of  the  saune  form  as 

Eqs  (103),  (104),  and  il06)  may  also  be  obtained  by  eliminating  the 
variable  Q  instead  of  p  in  Eqs  (85)  and  (86).  The  procedure  is  exactly 
the  same  as  previously  shown.  These  complementaiy  equations  are 
summarized  in  Appendix  E. 
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VI.  Nujoerical  Results 


Summary  of  Figures 


A  summary  of  the  performance  of  the  finite-difference  model 
developed  in  Chapters  II  through  IV  is  shown  in  Figures  14  through  70. 
The  capabilities  and  limitations  of  this  model  are  demonstrated  for  the 
four  basic  transmission  line  configurations  listed  in  Table  5. 


Table  5 

TRANSMISSION  LINE  CONFIGURATIONS  ANALYZED 
WITH  FINITE-DIFFERENCE  MODEL 


LINE 

CONFIGURATION 

ID  (in) 

LENGTH  (ft) 

TERMINATION 

1 

0.305 

40 

open 

2 

0. 170 

5 

blocked 

3 

0. 170 

5 

volume 

4 

0. 170 

5 

open 

The  model  is  designed  to  predict  the  transient  response  of  fluid 
transmission  lines  for  input  signals  of  arbitrary  shape.  However,  only 
three  input  signals,  a  trapezoidal  pulse,  a  sinusoidal  pulse,  and  a 
terminated  ramp,  were  used  to  study  model  performance.  Input  signal 
characteristics  are  shown  in  Table  6.  The  pulses  provide  sudden 
acceleration  followed  by  sudden  deceleration  of  the  fluid  with  the 
trapezoidal  pulse  posing  the  most  difficulty  for  the  algorithm.  The 
terminated  ramp  input  was  used  in  lieu  of  a  step  function  since  the 


Tab3e  6 


Input  Signal  auid  Line  Characteristics  for  Figures  14-52 


Fluid  Properties 

Line  Characteristics 

lb 

H  »  0.044  - 

m 

L  = 

40  ft 

ft  hr 

D  = 

0 . 305  1 n 

lb 

M  = 

4 

p  =  0.0747 

m 

3 

r  = 

0.12  In 

ft 

1 

r  = 

0.14  in 

2 

T  =  530°  R 

r  = 

0.15  In 

0 

3 

r  - 

4 

0.1525  In 

c  =  1128.4 

ft 

Gridpoints  see  text 

0 

sac 

SIGNAL  CHARACTERISTICS 

SIGNAL 

RISE  TIME  DURATION  DECAY  TIME 

Trapezoidal 

0.0002  sec  0.005 

sac  0.0002  sec 

Pulse 

Sinusoidal 

-  0.005 

sec  — 

Pulse 

(X»5.64  ft) 

Step 

0.0002  sec  - 

—  - 

finite-difference  algorithm  does  not  cope  well  with  discontinuities. 
These  three  signals  provide  a  good  assessment  of  the  capabilities  of 
the  algorithm.  Additionally,  some  published  data  for  the  step  response 
of  fluid  transmission  lines  is  available.  'Comparisons  with  published 
data  are  made  where  possible. 

A  linearized  version  of  the  finite-difference  algorithm  was  also 
run  for  comparison  with  the  nonlinear  algorithm.  Comparison  of  linear 
amd  nonlinear  results  from  the  same  basic  algorithm  highlights  the 
capability  of  the  model  to  identify  nonlinear  effects  resulting  from  an 
input  signal  of  finite  amplitude. 
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In  addition  to  the  finite-difference  model  results.  Figures  71 


through  86  demonstrate  the  performance  of  the  modal  superposition 
method  developed  in  Chapter  V.  Comparisons  between  modal  superposition 
results  and  results  from  the  finite-difference  algorithms  are  also 
shown.  Figures  shown  in  this  chapter  are  summarized  in  Table  7.  The 
term  steady  friction  is  used  to  denote  friction  losses  calculated  using 
a  constant  Darcy-Weisbach  friction.  Unsteady  friction  denotes  viscous 
losses  calculated  directly  from  the  transient  velocity  profile  as 
opposed  to  a  constant  friction  factor.  In  addition  to  unsteady 
friction,  there  is  unsteady  heat  transfer  which  is  calculated  from  the 
transient  temperature  profile.  This  nomenclature  will  be  used  to 
simplify  discussion  throughout  this  chapter. 


Simulation 


Linear  and  nonlinear  versions  of  the  nximerical  algorithm  were  run 
using  1  psi  amplitude  input  signals  applied  at  x  =  0  in  the 
computational  domain  simulating  a  semi-infinite,  0.305  in  ID,  fluid 
transmission  line.  This  is  the  first  configuration  listed  in  Table  5. 
The  actual  length  of  the  computational  domain  simulated  for  this 
transmission  line  is  40  ft.  The  traveling  wave  produced  by  the  Input 
pressure  pulse  never  reaches  the  opposite  end  of  the  computational 
domain  during  the  time  period  of  interest.  This  prevents  any  portion 
of  the  traveling  wave  from  being  reflected,  producing  results 
approximating  the  transient  response  of  a  semi-infinite  transmission 
line.  The  effects  of  numerical  dispersion  and  dissipation  produced  by 
the  algorithm  as  well  as  the  effects  of  the  viscous  damping  and  heat 
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Table  7 


SUMMARY  OF  FIGURES 
Input 

FlgureLlneAlgorlthffi  Signal  Damping  Comments 


14 

0 

N 

TP 

None 

Profile,  40,  120,  and  200  Grid 

15 

0 

N 

TP 

None 

Time  History,  40,  120,  200  Grid 

16 

0 

L 

TP 

None 

Profile,  40,  120,  and  200  Grid 

17 

0 

L 

TP 

None 

Time  History,  40,  120,  200  Grid 

18 

0 

N/L 

TP 

None 

Time  History,  Fig  6.2,  Fig  6.3 

19 

0 

N 

SP 

None 

Profile,  40,  120,  and  200  Grid 

20 

0 

N 

SP 

None 

Time  History,  40,  120,  200  Grid 

21 

0 

L 

SP 

None 

Profile,  40,  120,  and  200  Grid 

22 

0 

L 

SP 

None 

Time  History,  40,  120,  200  Grid 

23 

0 

N 

TP 

None 

Profile/Time  History,  Filtering 

24 

0 

N 

TP 

SF 

Profile/Time  History 

25 

0 

L 

TP 

SF 

Profile/Time  History 

26 

0 

N 

TP 

UF 

Profile/Time  History 

27 

0 

L 

TP 

UF 

Profile/Time  History 

28 

0 

N/L 

TP 

UF 

Profile/Time  History 

29 

0 

N 

TP 

SF/UF 

Profile/Time  History 

30 

0 

L 

TP 

SF/UF 

Profile/Time  History 

31 

0 

N 

TP 

FH 

Profile/Time  History 

32 

0 

L 

TP 

FH 

Profile/Time  History 

33 

0 

N/L 

TP 

FH 

Profile/Time  History 

34 

0 

N 

TP 

UF/FH 

Profile/Time  History 

35 

0 

L 

TP 

UF/FH 

Profile/Time  History 

36 

0 

N 

TP 

SF/FH 

Profile/Time  History 

37 

0 

N 

TR 

UF 

Profile/Time  History 

38 

0 

L 

TR 

UF 

Profile/Time  History 

39 

0 

N 

TR 

FH 

Profile/Time  History 

40 

0 

L 

TR 

FH 

Profile/Time  History 

41 

0 

N/L 

TR 

UF&FH 

Profile 

42 

0 

N/L 

TR 

UF&FH 

Tlri’e  History 

43 

0 

N&L 

TR 

UF/FH 

Profile 

44 

0 

N&L 

TR 

UF/FH 

Time  History 

A5 

0 

N 

SP 

FH 

Profile,  1,  5,  10  psi  amplitude 
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Table  7  (continued) 


SUMMARY  OF  FIGURES 


Ramp 


Input 


Figure  Line  Algorithm  Signal  Damping 


Comments 


46 

0 

L 

SP 

FH 

47 

0 

N/L 

SP 

FH 

48 

0 

N/L 

SP 

FH 

49 

0 

N/L 

SP 

FH 

50 

0 

L 

TP 

FH 

51 

0 

L 

TP 

FH 

52 

0 

L 

TR 

ni 

53 

- 

- 

— 

— 

54 

B 

L 

TR 

FH 

55 

B 

L 

TR 

FH 

56 

B 

N/L 

TR 

FH 

57 

V 

L 

TR 

FH 

58 

V 

L 

TR 

FH 

59 

V 

L 

TR 

FH 

60 

V 

L 

TR 

FH 

61 

V 

L 

TR 

FH 

62 

0 

N/L 

TR 

FH 

63 

0 

L 

TR 

FH 

64 

0 

L 

TR 

FH 

65 

0 

L 

TR 

FH 

iS 

w’ 

i_ 

TR 

FH 

67 

3/0 

L 

TR 

FH 

68 

B/0 

L 

TR 

FH 

69 

- 

- 

— 

— 

70 

SC 

L 

— 

Ur 

71 

B 

L 

TR 

— 

72 

B 

L 

TR 

— 

73 

B 

L 

TR 

— 

74 

B 

L 

TR 

— 

75 

B 

L 

TR 

— 

78 

B 

L 

TR 

— 

77 

B 

L 

TR 

— 

Profile,  1,  5,  10  psi  amplitude 
Profile,  1,  5,  10  psi  amplitude 
Time  History,  1,  5,  10  psi 
Convergence  of  N  and  L  Results 
Laminar/Turbulent  Damping 
Different  Interval  Boundaries 
Velocity  Profile  Snapshots 
Line  Configurations 
Laminar/Turbulent  Friction,  Dat? 
Time  History,  0.25,  1.0,  2.0  psi 
Time  History 

Time  History,  0.25,  1.0,  2.0 

Time  History,  0.5  V^,  Inlet 

Time  History,  V^=  0.5  V^,  Exit 

Time  History,  V  *  1.0  V  ,  Inlet 
c  ^ 

Time  History,  V^=  1.0  V  ,  EIxit 
Time  History,  0.25  psi 
Time  History,  1.0  psi.  Inlet 
Time  History,  1.0  psi.  Exit 
Time  History,  2.0  psi.  Inlet 
T’me  1.0  psi.  Exit 

Time  History,  i.O  psi.  Inlet 
Time  History,  2.0  psi.  Inlet 
Polynomial  Velocity  Profiles 
Oil-Filltd  Line,  rinsure 

Steady  LRC/Finite-Difference 
Unsteady  LRC/Finite-Difference 
Unsteady  LRCG/Finite-Difference 
Unsteady  LRCGAJnsteady  LRC 
Unsteady  LRCG,  6  Modes 
Unsteady  LRCG,  10  Modes 
Unsteady  LRCG,  15  Modes 


Tab 1 e  7  ( cont 1 nued ) 


SUMMARY  OF  FIGURES 
Input 


Figure 

Line 

Algorithm 

Signal  Damping 

Comments 

78 

B 

L 

TR 

— 

Unsteady  LRCG,  20  Modes 

79 

V 

L 

TR 

— 

LRCG/Finite-Difference,  Inlet 
15  Modes,  V  =  V  ,  Pressure 

C  L 

80 

V 

L 

TR 

— 

LRCG/Finite-Difference,  Inlet 
15  Modes,  Velocity 

81 

V 

L 

TR 

— - 

LRCG/Finite-Difference,  Exit 
15  Modes,  V^  “ 

82 

B 

L 

TR 

— 

LRCG,  Branched  Line 

83 

0 

L 

TR 

— 

LRCG/Flnite  Difference,  Inlet 
15  Modes,  Velocity,  1.0  psi 

84 

0 

L 

TR 

LRCG/Flnite  Difference.  Exit, 
15  Modes,  Velocity,  1.0  psi 

85 

0 

L 

TR 

— 

LRCG/Finite  Difference,  inlet 
15  Modes,  Velocity,  0.25  psi 

86 

0 

L 

TR 

— — 

LRCG/Flnite  Difference,  Exit, 
15  Modes,  Velocity,  0.25  psi 

B  = 

Blocked  Line 

TP  = 

Trapezoidal  SF  =  Steady  Friction 

0  * 

Open  Line 

SP  = 

Sinusoidal  UF  *  Unsteady  Friction 

V  = 

Volume 

Terminated 

TR  = 

Terminated  FH  =  Combined  Unsteady 

Line  Ramp  Friction  and  Heat 

Transfer 


transfer  functions  can  be  more  easily  observed  with  no  reflections 
present. 

Figure  14  shows  the  profiles  of  a  trapezoidal  pressure  pulse  at 
three  separate  times  as  it  travels  from  left  to  right  along  a  line  with 
no  damping.  The  nonlinear  version  of  the  numerical  algorithm  for  Eqs 
(54),  (55),  (60),  and  (61)  is  used  with  the  viscous  damping  and  heat 
transfer  terms  neglected.  As  discussed  in  Chapter  II,  the  damped-wave 
equations  are  coupled  together  through  the  damping  and  heat  transfer 
functions.  Setting  these  functions  to  zero,  as  in  Figure  14,  is 
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equivalent  to  specifying  a  uniform  velocity  over  the  cross  section, 
making  the  system  of  equations  equivalent  to  a  one-dimensional  system. 
In  the  nonlinear  algorithm,  coefficients  containing  axial  velocity,  u, 
amd  sonic  .elocity,  c,  vary  locally,  with  c  being  a  function  of  the 
local  temperature  in  this  model.  In  contrast,  the  linear  version  of 
the  algorithm  neglects  the  axial  velocity  in  these  coefficient  terms 
and  uses  a  constant  value  of  c  based  on  the  line  conditions  present 
before  signal  application.  The  three  curves  shown  in  each  part  of 
Figure  14  represent  a  plot  of  pressure  versus  aocial  distance  at  three 
separate  times.  These  are  three  snapshots  of  the  same  pressure  pulse 
as  it  travels  down  the  transmission  line.  This  figure  highlights 
effects  caused  by  numerical  errors  introduced  by  the  algorithm  and 
shows  the  importance  of  choosing  a  sufficiently  small  value  for  the 
spatial  discretization. 

Figure  14a  was  produced  using  a  computational  domain  containing  40 
grid  points  over  a  distance  of  40  ft  (1  grid  point  per  ft).  This 
severe  rounding  results  from  the  use  of  an  insufficient  number  of 
points  to  accurately  convey  the  information  contained  in  the  input 
signal  to  the  algorithm.  The  input  signal  that  the  algorithm  "sees" 
lacks  most  of  the  detail  of  the  original  input  signal.  The  high 
frequency  components  necessary  to  resolve  the  corners  and  steep  slopes 
of  the  pulse  cannot  be  carried  by  the  coarse  mesh. 

In  addition  to  the  distortion  caused  by  a  lack  of  waveform  detail 
for  the  input  signal,  numerical  distortion  appears  in  the  form  of 
oscillations  at  the  leading  edge  of  the  waveform.  These  oscillations 
result  from  the  niimerical  dispersion  error  introduced  by  the  algorithm. 
The  oscillations  are  not  very  severe  in  Figure  14a  since  Ax  is  quite 
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large  for  the  coarse  grid.  The  value  of  the  error  term  containing  the 
fourth  derivative  in  the  modified  equation  Eq  (80)  becomes  significant, 
providing  numerical  dissipation.  This  tends  to  damp  out  the 
oscillations  produced  by  the  third  derivative  error  term. 

Figures  14b  auid  c  use  120  and  200  grid  points,  respectively, 
providing  significantly  Improved  resolution  of  the  waveform.  The 
oscillations  produced  by  the  numerical  dispersion  become  more  severe 
due  to  a  decrease  in  the  damping  provided  by  the  fourth  derivative 
error  term  of  the  modified  equation.  The  oscillations  in  Figure  14c 
are  very  severe,  and  actually  extend  well  below  the  horizontal  axis  of 
the  chart.  The  graphics  program  used  to  produce  the  chart  ignores 
points  with  pressure  values  less  than  the  lower  limit  of  the  vertical 
axis,  exaggerating  the  distortion  to  the  leading  edge  of  the  waveform. 

Notice  that  the  slope  of  the  leading  edge  of  the  waveform  steepens 
while  the  slope  of  the  trailing  edge  flattens  as  it  progresses  down  the 
line.  This  is  a  nonlinear  effect  caused  by  locally  varying  values  of  c 
and  u  in  the  coefficient  terms  of  the  damped  wave  equations.  Results 
from  the  linear,  finite-difference  algorithm  show  no  steepening  of  the 
waveform,,  as  will  be  seen  in  Figures  16  and  17. 

Figure  15  shows  the  signal  pressure  as  a  function  of  dimensionless 
time  at  the  0,  10,  and  20  ft  tube  axial  locations.  This  is  the  same 
trauismission  line  shown  in  Figure  14  but  the  waveform  information  is 
shown  at  discrete  locations  as  a  function  of  time  instead  of  discrete 
times  as  a  function  of  location.  The  same  waveform  shape  change  and 
oscillations  observed  in  Figure  14  appear  here.  Note  that  the  curves 
in  Figure  15a  appear  much  cruder  than  their  counterparts  in  Figure  14a. 
In  Figure  14a,  the  pressure  is  mapped  at  each  grid  point  contained  in 
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r  ~ 


8.  40  Gi'id  Points 


b.  120  Grid  Points 


c.  200  Grid  Points 


Fig.  14.  Single  Trapezoidal  Pulse  Traveling  to  Right,  Nonlinear 
Results,  No  Damping,  Times  and  Relative  Positions  of  Leading  and 
Trailing  Edges  Are  The  Same  for  All  Signal  Profiles. 
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the  30  ft  span  shovm  on  the  horizontal  axis.  In  Figure  15a,  the 

Av 

pressure  is  calculated  at  each  value  of  time,  with  At  =  —  For 
long  lines  with  relatively  few  grid  points.  At  may  be  quite  large 
compared  to  the  rise  and  decay  times  of  the  waveform.  With  more  grid 
points,  the  resolution  increases  dramatically  as  shown  in  Figures  15b 
and  15c.  These  results  clearly  show  that  the  choice  of  computational 
grid  resolution  must  consider  the  waveform  of  the  input  signal.  Slowly 
changing  signals  require  fewer  grid  points  than  signals  with  steep 
slopes  amd  spikes  to  achieve  good  resolution.  Note  that  there  is  no 
significant  change  in  the  amplitude  of  the  wave  in  Figures  14b,  14c, 
ISb,  and  15c.  Some  numerical  dissipation  will  be  present  using  the 
nonlinear  form  of  the  algorithm  since  the  local  Courant  number  v  is  not 
always  exactly  1.0,  however,  it  causes  no  noticeable  damping  effect  on 
the  waveform  over  the  time  period  of  interest.  When  v  =  1  ,  the  error 
terms  in  the  modified  equation  Eq  (78)  for  the  numerical  algorithm  are 
zero.  No  oscillations  at  the  leading  edge  of  the  wave  are  present  when 
the  linear  version  of  the  algorithm  is  used.  The  value  of  v  is  1.0  at 
all  points  and  all  times  for  the  linear  algorithm.  For  a  fixed, 
uniform  grid,  the  value  of  v  cannot  always  be  1.0  with  the  nonlinear 
algorithm.  As  a  result,  numerical  errors  appearing  as  oscillations  in 
the  solution  are  produced. 

Results  from  application  of  the  numerical  algorithm  to  Eqs  (10) 
through  (12),  the  linear,  damped-wave  equations,  are  shown  in  Figures 
16  and  17.  The  same  semi-infinite  line  used  in  Figures  14  and  15  is 
used  here.  The  effect  of  the  number  of  grid  points  on  the  resolution 
can  again  be  seen,  however  the  waveform  undergoes  no  shape  change  as  it 
travels  down  the  tube.  Also,  for  the  reasons  discussed  above,  the 
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TIME  <tCo/g 

a.  40  Grid  Points 


TIME  <tCo/g 

b.  120  Grid  Points 


TIME  <tCo/g 

c.  200  Grid  Points 

Fig.  15.  Time  History  of  Single  Trapezoidal  Pulse  at  0,  10,  20  ft. 
Nonlinear  Results,  No  Daoiping,  Relative  Position  of  Leading  and 
Trailing  Edges  Are  The  Same  for  All  Signal  Time  Histories. 
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DISTANCE  (FEED 
a.  40  Grid  Points 


b.  120  Grid  Points 


c.  200  Grid  Points 


Fig.  16.  Single  Trapezoidal  Pulse  Traveling  to  Right,  Linear  Results, 
No  Damping 
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PRESSURE  (PSIA)  PRESSURE  (PSIA)  PRESSURE  (PSIA) 


a.  40  Grid  Points 


TIME  <tCo/L) 

b.  120  Grid  Points 


c.  200  Grid  Points 


Fig.  17.  Time  History  of  Single  Trapezoidal  Pulse  at  0,  10,  20  ft 
Linear  Results,  No  Damping 


oscillations  are  absent  at  the  leading  edge  of  the  waveform.  The  loss 
of  high-frequency  Information  with  a  coarse  grid  is  still  evident. 
Figure  18  directly  compares  the  linear  auid  nonlinear  results  from  the 
previous  figures.  Although  the  input  signal  amplitude  of  1  psi  is  less 
thaui  seven  percent  of  the  Initial  14.67  psia  undisturbed  condition 
there  is  noticeable  difference  between  the  linear  and  the  nonlinear 
results.  As  the  input  signal  amplitude  decreases,  the  linear  and 
nonlinear  results  will  converge.  Conversely,  higher  input  signal 
aunplltudes  will  cause  the  nonlinear  results  to  deviate  significantly 
from  linear  results.  This  will  be  shown  in  Figures  47  through  49. 

A  1  psi  sinusoidal  pulse  is  used  in  Fig^lres  19  through  22.  The 
input  pulse  is  the  positive  half  of  a  sinusoid.  The  pulse  begins  at  an 
amplitude  of  0  and  peaks  at  1  psi.  With  a  shallower  slope  than  the 
trapezoidal  pulse,  the  distortion  present  at  the  leading  edge  of  the 
waveform  for  the  nonlinear  case  shown  in  Figures  19  and  20  is  not 
nearly  as  severe.  However,  the  shape  change  of  the  waveform  as  it 
travels  down  the  tube  is  obvious  with  the  oscillations  appearing  at  the 
leading  edge  as  the  wave  form  steepens.  This  waveform  evolution  is 
consistent  with  experimental  results  obtained  by  Sandoz  (59).  The 
roughly  sinusoidal  pulses  used  by  Sandoz  exhibit  the  same  steepening  of 
the  leading  edge  of  the  waveform  shown  by  the  numerical  algorithm 
developed  in  this  study.  Figures  21  auid  22  show  the  linear  results  for 
this  case. 

The  oscillations  at  the  leading  edge  of  the  waveform  may  be 
filtered  out  by  a  simple  algorithm  used  to  post-process  the  results  of 
the  primary  finite-difference  algorithm.  An  example  of  a  simple 
smoothing  algorithm  described  by  Vliegenthart  is  shown  below  (72): 


93 


a.  40  C 


b-  120 


c.  200  GrJ 


Fig.  19.  Single  Positive  Half-! 
Nonlinear  Results,  No  Damping 
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PRESSURE  (  PSIA)  PRESSURE  (  P8IA)  PRESaRE  (  P 


IB.  20 


14.00 


14.40 


14,00^ -f- 
0.0 


Linear  Results,  No  Damping 


TIME  <tCo/L) 

a.  40  Grid  Points 


TIME  <tCo/g 

b.  120  Grid  Points 


TIME  <tCo/g 

c,  200  Grid  Points 

Fig.  22.  Time  History  of  Single  Sinusoidal  Pulse  at  0,  10,  20  ft. , 
Linear  Results,  No  Damping 
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(130) 


'  1 

q  =  T  +  2q  +  q  ) 

4  ^1^1  i-i 

I 

q^  is  replaced  by  q^  at  each  tlmestep.  If  the  post-processing 
algorithm  is  not  chosen  carefully,  attenuation  of  the  higher  frequency 
components  present  in  the  solution  may  occur,  eliminating  rapidly 
fluctuating  portions  of  the  solution.  Equation  (130)  maintains  formal 
second  order  accuracy,  but,  as  can  be  seen  from  Figure  23,  introduces  a 
significant  numerical  error  which  appears  as  rounding  of  the  sharp 
corners  of  the  original  trapezoidal  waveform.  The  loss  of  important 
waveform  features  appears  far  more  serious  than  auiy  error  resulting 
from  the  oscillations.  Other  more  sophisticated,  higher-order, 
filtering  algorithms  are  available:  however,  they  require  a  significant 
increase  in  computer  resources  to  implement.  They  must  be  applied  at 
every  grid  point  although  not  necessarily  at  each  time  step  (16:13; 

38: 168): 


'  . 
q  =  q  +  1/  __  A 

1  Ax 


Aq  Aq 


(131  ) 


Equation  (131)  is  third-order  accurate  and  does  reduce  numerically 
induced  oscillations  without  significantly  attenuating  important 
information  (38:168).  However,  when  performed  at  each  grid  point  and 
every  few  time  steps,  it  significantly  increases  the  number  of 
arithmetic  wpe'-ations  required  to  obtain  the  solution.  As  can  be  seen, 
use  of  filtering  algorithms  may  introduce  error  into  the  solution  as 
well  as  increase  the  computational  workload  for  the  computer.  The 
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14.  40 


14.20  - 
^4,00- 


0.0  6.0 


10.0  15.0  20.0  25.0  30.0 

distance  (FEED 


a.  Profiles  at  0.0046,  0.014,  and  0.023  sec 


b.  Time  History  at  0,  10,  and  20  ft 


Fig.  23.  Single  Trapezoidal  Pulse,  Nonlinear  Results,  No  Damping, 
Filtering  Algorithm  Applied,  200  Gi  iu 
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finite-difference  model  developed  In  the  present  study  is  part  of  a 
design  tool  meant  for  use  on  modest  computers  easily  accessible  to 
designers  of  fluid  systems.  With  this  goal  in  mind,  computer  resources 
required  by  this  model  are  being  kept  as  small  as  possible.  Use  of 
this  model  on  modern  supercomputers  would,  of  course,  allow  better 
filtering  algorithms  and  other  difference  operators  to  be  used.  No 
higher-order  filtering  algorithms  are  used  or  recommended  for  this 
particular  finite-difference  model  in  order  to  limit  computer 
processing  requirements. 

Viscous  damping  using  the  Darcy-Weisbach  friction  factor  has  been 
added  to  the  nonlinear  algorithm  for  Figures  24  and  25.  The  algorithm 
Is  still  equivalent  to  a  one -dimensional  algorithm,  as  in  the  previous 
figures,  since  the  damping  is  calculated  as  a  function  of  average 
velocity.  It  is  important  to  note  that  the  steady-state  friction 
results  shown  in  Figures  24  and  25  were  obtained  using  a  constant 
friction  factor  f  based  on  the  predicted  peak  value  of  transient 
average  velocity  in  the  transmission  line  with  no  damping  present. 

This  method  was  used  for  the  present  study,  since  calculating  f  as  a 
function  of  the  transient  average  velocity  produced  no  discernible 
attenuation  in  the  numerical  result.  This  was  also  the  approach  used 
by  Sandoz  (59).  The  attenuation  produced  by  steady-state  friction  will 
be  compared  to  the  attenuation  produced  by  unsteady  friction  in  Figures 
29  and  30. 

The  next  several  figures  employ  the  damping  expression  shown  in  Eq 
(16),  modified  for  transient  flow.  u(th  the  addition  of  the  unsteady 
friction  term,  velocity  profiles  may  be  obtained,  producing  a  solution 
containing  a  two-dimensional  velocity  field.  Recall  that  the  cross 
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section  of  the  transmission  line  is  broken  into  four  intervals  where  an 
average  velocity  is  calculated  for  each  interval.  The  intervals  are 
narrow  near  the  wall  to  provide  a  higher  resolution  where  the  velocity 
changes  are  the  most  extreme.  The  intervals  widen  toward  the  center  of 
the  transmission  line. 

The  velocity  profile  is  recovered  from  these  average  velocities 
and  used  to  calculate  the  velocity  gradients  for  use  with  Eq  (16). 
Figures  26  and  27  show  results  from  the  nonlinear  and  linear  algorithms 
with  unsteady  frictional  damping  for  a  1  psi  trapezoidal  pulse  input 
signal.  In  both  the  linear  and  the  nonlinear  cases  there  is  a 
significant  rounding  of  the  leading  edge  of  the  pulse.  Figures  26  and 
27  show  that  the  unsteady,  frictional  damping  produces  a  significantly 
different  waveform  than  steady-state  frictional  damping.  This  is 
supported  by  the  results  of  Wood  and  Funk  (74).  Also,  for  the 
nonlinear  case,  the  effect  of  unsteady  damping  on  the  waveform  appears 
to  overwhelm  numerical  distortion  present  at  points  away  from  the 
leading  edge  based  on  comparisons  with  Figures  14  and  15  where  no 
damping  is  present.  The  effect  of  the  oscillations  on  most  of  the 
waveform  does  not  appear  to  be  significant  in  the  presence  of  the 
unsteady  damping  at  this  input  sigi^ai  level  and  the  propagation 
distances  considered. 

Results  from  both  linear  and  nonlinear  algorithms  with  unsteady 
friction  are  directly  compared  in  Figure  28.  Differences  in  waveform 
result  from  shape  evolution  produced  by  the  nonlinear  algorithm.  Aside 
from  the  presence  of  the  numerical  oscillations  in  the  nonlinear 
results,  the  waveforms  are  very  similar  in  shape  and  magnitude. 

A  comparison  of  steady  versus  unsteady  frictional  damping  is  shown 
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Fig.  26.  Single  Trapezoidal  Pulse,  Nonlinear  Results,  Unsteady 
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Fig.  27.  Single  Trapezoidal  Pulse,  Linear  Results,  Unsteady 
Friction 
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Fig.  28.  Single  Trapezoidal  Pulse,  Comparison  of  Linear  (Fig.  27) 
and  Nonlinear  (Fig.  26)  Results 
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for  the  nonlinear  case  In  Figure  29  and  the  linear  case  in  Figure  30. 
Here,  both  steady  auid  unsteady  damping  functions  were  calculated  for 
the  sauae  0.305  in  ID  transmission  line.  The  attenuation  of  the 
waveform  produced  using  steady  friction  depends  on  the  choice  of  the 
friction  factor.  It  is  possible  to  select  a  value  for  the  friction 
factor  that  will  produce  a  desired  level  of  attenuation.  However,  this 
must  be  done  on  a  trial  and  error  basis  as  shown  by  Sandoz  (59).  The 
unsteady  friction  models  used  in  this  study  are  not  empirical  in  nature 
auid  do  not  require  a  trial  and  error  process  to  determine  accurate 
values  of  damping  as.  The  friction  factor  is  empirical  in  nature  and 
only  valid  for  flows  that  are  steady  or  very  near  steady  with  fully 
developed  velocity  profiles.  It  loses  its  validity  in  the  transient 
flows  induced  by  the  passage  of  a  pressure  wave. 

Figures  31  and  32  show  the  waveforms  resulting  from  the  nonlinear 
and  linear  algorithms,  respectively,  using  combined  unsteady  friction 
and  heat  transfer.  The  turbulent  component  of  the  total  thermal 
conductivity  is  calculated  using  Reynolds  auialogy  to  obtain  a  value 
for  k  from  u  .  The  laminar  component  of  k  is  the  physical 
property  k  of  the  fluid  under  consideration.  Figures  31  and  32  are 
combined  in  Figure  33,  again  showing  a  noticeable  difference  between 
the  linear  and  nonlinear  results.  Here,  the  signal  attenuation  for  the 
linear  results  is  greater  than  that  for  the  nonlinear  results.  In 
contrast.  Figure  28,  which  compares  linear  and  nonlinear  results  for 
unsteady  friction  alone,  shows  almost  no  difference  in  magnitude 
between  the  linear  and  the  nonlinear  results.  The  presence  of 
convective  terms  in  the  nonlinear  energy  equation  is  thought  by  the 
author  to  be  responsible  for  this  behavior.  The  only  difference 
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Fig.  30.  Single  Trapezoidal  Pulse,  Comparison  of  Steady  (Fig.  25) 
and  Unsteady  (Fig.  27)  Friction,  Linear 
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Fig.  31.  Single  Trapezoidal  Pulse,  Nonlinear  Results,  Combined 
Unsteady  Friction  and  Heat  Transfer 
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Fig.  33.  Single  Trapezoidal  Pulse,  Comparison  of  Linear  (Fig.  32) 
and  Nonlinear  (Fig.  31)  Results,  Combined  Unsteady  Friction  and  Heat 
Transfer 
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between  the  linear  and  nonlinear  versions  of  the  damped-wave  equations 
In  split  coefficient  form  are  the  factors  u  +  c  amd  u  -  c  appearing  in 


Eqs  (54)  and  (55).  In  the  linear  equations,  u  is  zero  and  c  is 

n 

constant.  For  the  signal  amplitudes  under  consideration  in  this  study, 
the  magnitude  of  u  +  c  and  u  -  c  in  Eqs  (54)  and  (55)  varies  only  a  few 

n  n 

percent  from  the  constant  value  c  used  in  the  linear  equations. 

However,  this  small  becomes  signiflcamt  for  signals  of  linger  duration 

as  will  be  seen  when  the  terminated  ramp  signal  is  considered.  The 

linear  and  nonlinear  energy  equations  used  to  calculate  the  temperature 

3T  2  a  f 1  1 

field  differ  with  the  presence  of  the  terms  u  -s-”  and  r  u  — 

n  3x  yc  n  3x 


P 

in  the  nonlinear  energy  equation  Eq  (7).  The  contribution  of  these 
terms  appears  significant  as  shown  by  Figure  33. 

A  direct  comparison  of  nonlinear  results,  with  and  without 
unsteady  heat  transfer,  for  the  trapezoidal  pulse  is  shown  in  Figure 
34.  The  waveform  obtained  using  combined  unsteady  friction  and  heat 
transfer  exhibits  more  rounding  at  the  leading  edge  amd  a  s’ightly 
higher  peak  at  the  trailing  edge  than  the  waveform  with  unsteady 
friction  alone.  The  peak  amplitudes,  however,  are  virtually  the  same. 

Figure  35  is  the  linear  analog  of  Figure  34  and  shows  a  definite 
difference  between  the  magnitudes  of  the  waveforms.  In  each  chart,  the 
curve  of  lesser  magnitude  results  from  combined  unsteady  friction  and 
heat  transfer  while  the  remaining  curve  is  produced  using  unsteady 
friction  only.  The  difference  between  waveforms  seen  in  Figure  34 
(more  rounding  at  the  leading  edge  and  a  slightly  higher  peak  at  the 
trailing  edge)  is  absent  for  the  linear  case  in  Figure  35.  This 
comparison  indicatas  the  importance  of  the  convective  term  in  the 
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Fig  34.  Single  Trapezoidal  Pulse,  Comparison  of  Unsteady  Friction 
(Fig.  26)  and  Combined  Unsteady  Friction  and  Heat  Transfer  (Fig.  31) 
Nonlinear  Results 
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energy  equation  with  regard  to  waveshape  and  attenuation.  A 
significant  transport  of  thermal  energy  by  convection  Is  apparent.  As 
a  point  of  reference,  nonlinear  steady  friction  and  unsteady  combined 
friction  and  heat  transfer  results  for  a  1  psl  trapezoidal  pulse  are 
compared  In  Figure  36. 

Pulses  have  been  examined  so  far,  but  the  step,  or  more  properly 
In  this  case,  the  terminated  ramp,  also  provides  Information  on  the 
behavior  of  these  algorithms.  Unlike  the  pulse,  the  terminated  ramp 
has  no  trailing  edge  to  produce  a  sudden  deceleration  of  the  fluid. 

The  effects  of  the  unsteady  turbulent  friction  model  Is  much  easier  to 
observe. 

Nonlinear  results  for  unsteady  frictional  damping,  alone,  and 
combined  friction  amd  heat  transfer  are  shown  in  Figures  37  and  39 
respectively.  Linear  results  are  shown  in  Figures  38  and  40.  The 
snapshots  at  successive  times  shown  in  these  four  charts  are  useful  in 
indicating  the  rate  of  daunping  as  the  wavefront  travels  along  the 
transmission  line.  A  nonlinear  effect  on  the  signal  amplitude  not  seen 
in  pulses  shows  up  with  the  terminated  ramp  input  signal  for  the  case 
of  unsteady  friction  alone.  Recall  that  for  the  pulse  in  Figure  28 
there  was  no  significant  difference  between  the  amplitude  of  the  linear 
and  nonlinear  results  for  the  trapezoidal  pulse  as  it  traveled  along 
the  tramsmission  line.  Only  the  shape  of  the  waveform  changed.  With 
the  terminated  ramp,  there  is  a  very  noticeable  difference  between  the 
linear  and  nonlinear  results  as  can  be  seen  in  Figure  41a  which 
compares  results  from  Figures  37a  and  38a.  The  curve  with  the  greater 
magnitude  in  each  pair  is  the  nonlinear  result.  This  difference 
between  waveforms  would  indicate  that  this  effect  increases  in 
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Fig.  36.  Single  Trapezoidal  PuJ 
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Fig.  37.  Terminated  Ramp,  Nonlinear  Results,  Unsteady  Friction 
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Fig.  39.  Terminated  Ramp,  Nonlinear  Results,  Combined  Unsteady 
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Fig.  40.  Terminated  Ramp,  Linear  Results,  Combined  Unsteady 
Friction  and  Heat  Transfer 


magnitude  as  pulse  duration  Increases.  This  difference  between  linear 
and  nonlinear  results  was  not  seen  In  the  trapezoidal  pulse  due  to  Its 
short  duration.  In  a  pulse,  the  flow  Is  cut  off  as  the  trailing  edge 
of  the  pulse  passes,  preventing  the  nonlinear  convective  terms  from 
contributing  significantly  to  the  amplitude.  This  would  suggest  that 
the  longer  the  pulse  duration,  the  greater  the  effect  of  the  nonlinear 
convection  terms  In  the  damped  wave  equations. 

Linear  and  nonlinear  results  with  combined  friction  and  heat 
transfer  are  compared  In  Figure  41b.  The  curve  with  the  greater 
amplitude  of  each  pair  Is  the  nonlinear  result.  This  same  relationship 
between  linear  amd  nonlinear  results  was  observed  in  Figure  33  for  the 
trapezoidal  pulse.  The  difference  between  the  linear  and  nonlinear 
results  Is  due  to  the  presence  of  the  nonlinear  terms  In  the  energy 
equation,  as  It  was  for  the  trapezoidal  pulse.  Figures  42a  and  42b  are 
the  time  histories  for  Figures  41a  amd  41b  respectively. 

Figures  43  and  44,  again,  show  the  presence  of  the  heat  transfer 
damping  term  produces  a  noticeable  difference  from  results  with 
frictional  damping  alone.  Linear  results  are  compared  in  part  a  of 
each  figure  while  nonlinear  results  are  compared  in  part  b.  In  each 
pair  of  curves,  the  one  with  the  greater  amplitude  represents  unsteady 
friction  alone.  The  remaining  curve  results  from  combined  friction  and 
heat  transfer. 

One  trend  absent  from  results  for  the  terminated  ramp  input  signal 
is  the  distortion  of  the  step  waveform  shown  experimentally  in  previous 
studies  of  long  transmission  lines  (5:7;  28:280;  30:500-501).  Braden 
amd  Kantola  noted  am  amplitude  dependent  effect  that  distorted  the 
leading  edge  of  the  step  at  points  downstream  of  the  point  of  signal 
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a.  Unsteady  Friction  (Figs  6.24a,  6.25a) 


b.  Combined  Unsteady  Friction  and  Heat  Transfer  (Figs  6.26a,  6.27a) 


Fig. 
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Terminated  Ramp,  Comparison  of  Linear  and  Nonlinear  Results 
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b.  Combined  Unsteady  Friction  and  Heat  Transfer  (Figs  6.26b,  6.27b) 


Fig.  42.  Terminated  Ramp,  Comparison  of  Linear  and  Nonlinear 
Results,  Time  History 
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b.  Nonlinear  Results  (Figs  6.25a,  6.27a) 


Fig.  43.  Terminated  Ramp,  Comparison  of  Unsteady  Friction  and 


Combined  Unsteady  Friction  and  Heat  Transfer 


application.  Kantola  described  this  effect  as  an  initial  rapid  rise 
followed  by  a  nearly  constant  amplitude,  persisting  for  some  time, 
followed  by  a  slower  rise.  He  attributed  the  inability  of  his  model  to 
predict  this  behavior  to  the  lack  of  nonlinear  terms  capable  of 
accounting  for  convected  mass,  momentum,  and  heat  transfer  in  his  basic 
equations.  Karam  and  Leonard  noted  similar  behavior  at  an  upstream 
location  in  a  50  ft  line  (30).  In  both  cases,  the  tubing  was  coiled. 
The  basic  equations  used  in  this  study  retain  the  nonlinear  convective 
terms  but  still  fail  to  predict  this  experimentally  observed  behavior. 
One  possible  explanation  is  that  nonlinear  convection  alone  is  not 
responsible  for  this  effect.  Additional  unsteady  frictional  effects 
caused  by  the  distortion  of  the  transient  velocity  profile  in  the 
coiled  tube  may  play  an  Important  role  in  producing  this  amplitude 
dependent  damping  effect.  No  data  for  comparison  on  straight  long 
lines  has  been  found  in  the  course  of  this  study.  Although  a  straight 
line  was  ",sed  by  Sandoz,  her  work  was  restricted  to  short  pulses  amd  no 
result  similar  to  that  described  above  was  rejxjrted. 

The  effect  of  increased  signal  amplitude  on  the  algorithm  is 
demonstrated  in  Figures  45  through  46  using  a  sine  pulse  as  am  input 
signal.  Figures  45  and  46  show  nonlinear  auid  linear  results 
respectively  for  combined  unsteady  friction  auid  heat  transfer.  More 
rapid  distortion  and  ultimate  formation  of  a  discontinuity  are  evident 
in  the  nonlinear  case  for  the  higher  pressures.  Over  the  30  ft 
distance  considered  in  the  graph,  the  higher  amplitude  input  signal 
formed  a  discontinuity  much  earlier  than  the  smaller  amplitude  inputs. 
Once  a  discontinuity  forms  using  this  algorithm,  the  accuracy  is 
suspect  since  this  algorithm  produces  severe  oscillations  which 
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Fig.  45.  Single  Positive  Half-Sinusoidal  Pulse  Traveling  to  Right, 
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Fig.  46.  Single  Sinusoidal  Pulse  Traveling  to  Right,  Linear  Results 
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overwhelm  the  rest  of  the  solution.  Also,  the  nonlinear  algorithm  Is 
invalid  in  the  presence  of  discontinuities  since  the  damped-wave 
eqxiations  are  in  nonconservative  form  and  may  permit  the  existence  of 
small  mass  sources  or  sinks  (2:52).  Special  treatment  of 
discontinuities  is  required.  When  the  linear  algorithm  is  used  the 
higher  input  signal  pressures  cause  no  oscillations  and  only  result  in 
signals  of  higher  amplitudes. 

Linear  and  nonlinear  results  are  compared  on  the  same  coordinate 
axes  in  Figures  47  and  48  for  sinusoidal  pulses  of  1  psi,  5  psl  and  10 
psl.  The  importance  of  nonlinear  effects  on  propagating  signals  is 
clearly  demonstrated  for  large  pressure  amplitudes. 

An  input  signal  amplitude  of  0.25  psi  produces  the  results  shown 
in  Figure  49.  The  behavior  of  the  nonlinear  algorithm  results  closely 
follow  those  produc€Ki  by  the  linear  algorithm.  Even  for  this  small 
signal  (1.7  percent  of  the  initial  pressure)  a  change  in  waveform  can 
still  be  observed  as  the  signal  travels  down  the  transmission  line. 

The  shape  distortion  is  a  cumulative  process  aind  small  effects 
accumulated  over  a  large  distance  become  significant. 

As  stated  before,  the  damping  functions  include  a  model  for 
turbulent  friction  which  engages  when  the  Reynolds  number,  based  on  the 
average  axial  velocity  at  each  axial  node  location,  exceeds  2300.  No 
data  indicating  criteria  for  the  onset  of  turbulence  in  transient  flow 
was  found.  The  value  of  2300  was  chosen  from  steady  flow  theory  for 
lack  of  a  better  value.  Figure  50  compares  laminar  and  turbulent 
results  for  the  same  transmission  line  configuration.  Linear  results 
were  used  for  this  comparison  so  that  nonlinear  effects  like  shape 
evolution  and  oscillations  vrould  not  mask  the  difference  between  the 
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Fig.  47.  Single  Positive  Half-Sinusoidal  Pulse  Traveling  to  Right, 
Comparison  of  Linear  (Fig.  45)  and  Nonlinear  (Fig.  46)  Results, 
Combined  Unsteady  Friction  and  Heat  Transfer 
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Fig.  48.  Single  Sinusoidal  Pulse,  Comparison  of  Linear  and  Nonlinear 
Results,  Combined  Unsteady  Friction  and  Heat  Transfer,  Time  History 
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Fig.  49.  Single  Positive  Half-Sinusoidal  Pulse,  Convergence  of 
Linear  and  Nonlinear  Results  for  Sufficiently  Small  Amplitude 
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caused  severe  instability  in  the  results.  This  problem  was  easily 
laminar  and  turbulent  friction  models.  The  waveform  with  the  greater 
amplitude  is  obtained  using  laminar  friction  only.  The  waveform 
produced  using  turbulent  friction  damping  has  a  steeper  slope  at  the 
leading  edge,  is  square  at  the  top,  and  decays  somewhat  faster  than 
the  laminar  result. 

So  far,  with  the  exception  of  Figures  24  and  25,  the  location  of 
the  boundaries  of  the  flow  Intervals  dividing  the  cross  section  of  the 
transmission  line  have  remained  fixed.  For  the  0.305  in  ID  line,  these 
boundaries  are  located  at  0.12,  0.14,  0.15,  and  0.1525  in  measured  from 
the  centerline.  Figure  51  compares  the  results  from  the  boundary 
locations  shown  above  to  results  from  boundaries  located  at  0.07625, 
0.122,  0.1375,  and  0.1525  in.  There  is  only  a  slight  difference  in  the 
waveform.  The  results  do,  however,  show  that  location  of  the  flow 
interval  boundaries  does  affect  the  results  to  a  small  degree.  The 
location  of  the  inner  boundary  of  the  outermost  flow  Interval  has  the 
most  effect  on  the  results.  Results  appear  to  be  insensitive  to  the 
locations  of  the  other  boundaries  if  the  thickness  of  the  outer 
interval  is  less  than  ten  percent  of  the  trsmsmission-line  radius. 

Large  fluctuations  in  the  value  of  the  damping  function  at  the 

wall  were  produced  when  too  narrow  an  outside  interval  (r  <  r  <  r  ) 

M*  1  M 

was  used.  These  fluctuations  often  fixed  by  decreasing  the  time  step, 
and  consequently  the  spatial  discretization,  with  an  associated 
increase  in  computer  cpu  time,  or  by  increasing  the  width  of  the 
outside  Interval.  To  ensure  good  resolution  at  the  wall,  the  inner 
boundary  of  the  outer  interval  was  kept  within  what  was  estimated  to  be 
the  viscous  sublayer  for  the  expected  flow  conditions.  This  type  of 
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a.  Profiles  at  0.0046,  0.014,  and  0.023  sec 


b.  Time  History  at  0,  10,  and  20  ft. 


Fig.  50.  Single  Trapezoidal  Pulse,  Comparison  of  Laminar  (Fig.  27) 


and  Turbulent  Unsteady  Friction 


b.  Tins  History  at  0,  10,  and  20  ft. 


Fig.  51.  Trapezoidal  Pulse,  Comparison  of  Two  Different  Sets  of  Flow 
Interval  Boundary  Locations 
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from  the  wall.  It  Is  the  effect  of  this  unsteady  velocity  profile  that 
Instability  has  also  been  noted  using  the  method  of  characteristics 
with  steady  damping  on  narrow  lines  (16:9).  More  detail  on  the  effects 
of  Interval  boundary  location  and  the  optimum  number  of  intervals  was 
not  pursued  in  this  study.  As  few  intervals  as  possible  were  used  in 
keeping  with  the  goal  to  minimize  the  computer  resources  required  to 
model  the  transmission  line. 

Use  of  the  piecewise  defined  polynomial  allows  information  about 
the  instantaneous  velocity  profiles  to  be  obtained  for  each  node 
location  at  each  time  step.  This  velocity  profile  Information  is  the 
basis  for  calculation  of  the  unsteady  friction  and  heat  transfer 
damping  functions  at  the  boundary  of  each  flow  interval.  Figures  52a 
through  52d  show  the  transient  velocity  profiles  in  a  line  at  the  point 
of  signal  application  and  at  a  point  5  ft  downstream.  Each  individual 
figure  is  a  snapshot  of  the  velocity  profile  caused  by  a  1  psl 
terminated  ramp  input  signal  as  it  p^'ogresses  down  the  transmission 
line.  Kct-  th#*  bulge  in  ths  profile  near  the  wall  as  the  flow  is 
initially  accelerated  by  the  arrival  of  the  pressure  wave  at  the  5  ft 
location.  Similar  behavior  of  the  velocity  profile  near  the  wall  was 
observed  experimentally  by  Shih  for  pressure  inputs  produced  by  means 
of  burst  discs  (63).  Significant  damping  occurs  near  the  wall  as  a 
result  of  this  bulge  in  the  velocity  profile  while  little  happens  away 
makes  the  use  of  the  friction  factor,  which  is  based  on  steady,  fully 
developed  flow,  inaccurate.  A  short  time  after  the  initial  arrival  of 
the  pressure  wave,  the  velocity  profile  begins  to  take  on  a  more 
familiar  shape  as  the  fluid  continues  to  be  accelerated  and  the  steady, 
velocity  profile  continues  to  develop. 
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a.  Flow  Accelerated  at  Station  1,  b.  Wavefront  Arrives  at  Station  2 
No  Signal  at  Station  2 


c.  Flow  Accelerating  at  Station  2  d.  Steady  Velocity  Prefile 

Developing  at  Station  2 

Fig.  52.  Instantaneous  Velocity  Profiles  Produced  by  Passage  of 
Terminated  Ramp  Input  Signal 


Lines  of  Finite  Length 


All  cases  discussed  so  far  Involve  approximations  of  semi-infinite 
transmission  lines.  This  configuration  kept  the  effects  of  numerical 
dispersion  and  the  damping  functions  on  the  Input  waveform  from  being 
masked  by  reflections  from  the  transmission  line  ends.  Finite  lines 
represent  more  realistic  configurations  often  found  in  complex 
pneumatic  and  hydraulic  systems.  An  important  goal  of  the  present 
study  is  to  arrive  at  a  practical  way  to  analyze  real  systems 
containing  one  or  more  transmission  lines. 

A  5  ft,  0.170  in  ID  air-filled  finite  tramsmission  line  was 
simulated.  Three  types  of  terminations  were  examined:  the  blocked 
line,  the  volume-terminated  line,  and  the  open  line.  These  are  shown 
in  Figure  53. 

Figure  54  shows  the  response  to  a  step  input  for  the  blocked  line 
configuration.  Results  with  turbulent  and  laminar  friction  are 
compared  to  data  obtained  by  Karam  auid  Leonard  (30:501).  The  curve 
with  the  most  severe  attenuation  was  produced  using  turbulent  friction. 
The  turbulent  and  launlnar  results  bracket  the  data  with  the  curve  shape 
following  the  shape  of  the  data  quite  well.  Figure  55  compares 
nonlinear  results  generated  by  terminated  ramp  inputs  of  0.25,  1.0,  and 
2.00  psi  for  the  5  ft,  0.170  in  ID  blocked  line.  The  magnitude  of 
peaks  of  the  transient  waveforms  differ  according  to  the  input  pressure 
amplitude  with  the  largest  input  pressure  producing  the  largest  peak. 
The  relation  among  the  amplitudes  of  the  peaks  is  retained  as  the 
signal  attenuates.  The  slope  of  leading  edge  of  each  waveform  steepens 
in  relation  to  the  input  pressure  amplitude,  with  the  steepest  slope 
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Fig.  53.  Finite-Length,  Transmission-Line  Configurations 


produced  by  the  largest  input  pressure.  Figure  56  compares  the  linear 

and  nonlinear,  finite-difference  results  for  a  2.  psl  input  pressure 

terminated  ramp.  Although  the  two  curves  are  in  close  agreement,  a 

definite  difference  in  slopes  can  be  seen.  A  5  ft  line  is  relatively 

short.  For  longer  lines,  the  difference  between  linear  and  nonlinear 

results  would  differ  by  considerably  more  due  to  the  longer  travel  time 

over  which  the  nonlinear  effects  could  accumulate. 

The  volume- terminated  line  is  examined  next.  Figure  57  compares 

the  pressure  response  to  a  step  input  for  the  volume-terminated  line 

with  V  equal  to  0.25V  ,  V  ,  and  2.0V  where  V  is  the  volume  of  the 
c  L  L  L  c 

terminating  chaunber  and  V^  is  the  volvime  of  the  5  ft  line.  Compression 
in  the  volume  is  assumed  to  occur  isothermal ly.  The  volume  is 
described  by  the  following  equation  (19:3): 


ap  jr  a  n  p  u 


(131) 


at 
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—  LAMINAR  FRICTION 


1.60 

1.40 

1.20 

1.00 

e.B0 

0.60 

0.40 

0.20 

0.00 


0.00  0.0t  0.02  0.03  0.04  0.05 


TIME  (SECONDS) 


Fig.  54.  Transient  Response  for  Blocked  5  ft,  0.170  in  ID  Line 
Linear  Results,  Terminated  Ramp  Input 


(p-p^)/(p. -p^) 


TIME  (SECONDS) 


Fig.  55.  Transient  Response  for  Blocked  5  ft.  0.170  in  TD  Line, 
Linear  Results,  Terminated  Ramp  Input 
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Fig.  56.  Comparison  of  Linear  and  Nonlinear  Response  to  2.0  psi 
Terminated  Ramp  in  5  ft,  0.170  in  ID  Blocked  Line 


PRESSURE  (PSIA) 


TIME  (tCo/L) 


Fig.  57.  1.0  psi  Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Volxane 


Terminated  Line 


where  n  =  1  (n  =  1.4  for  adiabatic  behavior).  The  volume  of  the 
termination  has  a  very  significant  effect  on  the  overall  response  of 
the  system.  The  behavior  of  the  average  line  inlet  and  exit  velocities 
is  shown  in  Figures  58,  59,  60  and  61.  Rapid  velocity  fluctuations 
with  sharp  reversals  can  be  seen. 

The  final  line  termination  examined  is  the  open-end  line  where  the 
end  of  the  tramsmission  line  is  maintained  at  ambient  pressure.  This 
configuration  proved  the  most  demanding  of  all  with  regards  to 
numerical  stability  of  the  algorithm.  A  comparison  of  linear  and 
nonlinear  results  is  shown  in  Figure  62  for  a  0.25  psi  terminated  ramp 
input.  The  results  are  identical  up  to  the  point  where  the  nonlinear 
results  diverge  and  continue  to  Increase  while  the  linear  results 
approach  a  steady-state  value.  Note  the  presence  of  periodic  velocity 
spikes  appearing  at  regular  Intervals  in  the  linear  results.  These 
spikes  are  more  pronounced  at  an  input  pressure  of  1.0  psi  as  shown  in 
the  linear  results  of  Figure  63  and  64.  These  spikes  are  thought  to  be 
the  main  cause  of  the  instability  of  the  nonlinear  numerical  algorithm 
when  used  with  the  open  line  configuration. 

Input  pressures  greater  than  0.25  psi  caused  the  nonlinear 
algorithm  to  fail  completely  unless  the  time  step,  and  consequently, 
the  computational  mesh  size  in  the  axial  direction  were  significantly 
reduced  in  size.  No  attempt  was  made  to  see  if  these  spikes  ultimately 
damp  out  or  remain  in  the  solution  indefinitely.  A  constant  pressure 
is  used  for  the  boundary  condition  at  the  open  end  of  the  line.  It  is 
still  unknown  whether  these  velocity  spikes  are  an  artifact  of  the 
numerical  algorithm  due  to  an  inaccurate  description  of  the  boundary 
condition  or  represent  a  phenomenon  that  could  be  experimentally 


145 
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Fig. 


62. 


Comparison  of  Linear  and  Nonlinear  Results  for  0.25  psi 


Terminated  Ramp  Input  to  5  ft,  0.170  in  ID  Open  Line 


50 


VELOCITY  AT  X  =  L  (FPS) 


observed.  No  experimental  or  theoretical  transient  velocity  data  for 
open-end  lines  was  found. 

The  linear  algorithm  was  used  to  examine  the  trauisient  response 
for  input  pressures  above  0.25  psi  since  the  nonlinear  algorithm  became 
unstable.  When  the  amplitude  of  the  input  pressure  step  was  increased 
to  2  psi,  the  velocity  spikes  became  almost  nonexistent  as  shown  in 
Figures  65  and  66.  Since  this  spike  phenomena  is  present  in  the  linear 
results  as  well  as  the  nonlinear  results,  it  cannot  be  solely 
attributed  to  the  nonlinear  algorithm. 

It  is  informative  to  look  at  a  comparison  of  the  inlet  velocity  of 
the  closed-end  and  open-end  line  configuration.  Figure  67  compares  the 
velocity  time  history  resulting  from  a  1.0  psi  input  pressure  using  the 
linear  algorithm.  Notice  that  up  until  t  =  2.0,  the  velocity  traces  of 
both  configurations  are  exactly  the  same.  The  same  is  true  in  Figure 
68  for  the  2  psi  input  pressure.  In  each  figure,  the  velocity  spikes 
and  steep  rises  occur  approximately  2.0  dimensionless  time  units  apart 
which  is  the  time  it  takes  the  wave  to  travel  twice  the  length  of  the  5 
ft  line.  The  spikes  are  obviously  produced  and  maintained  by  the 
downstream  open-end  boundary  condition.  It  should  also  be  noted  that 
the  respKsnse  of  the  transmission  line  is  exactly  the  same  until 
sufficient  time  passes  for  the  effects  of  the  downstream  boundary 
condition  to  propagate  to  the  inlet. 

All  cases  exaunlned  so  far  have  been  for  a  line  with  no  steady  mean 
flow.  One  advaintage  of  using  a  numerical  approach  to  model  a 
transmission  line  is  the  ease  of  incorporating  an  established  mean  flow 
as  a  set  of  initial  conditions.  However,  it  is  important  that  the 
method  used  to  obtain  the  initial  flow  field  lb  consistent  with  the 
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Fig.  67.  Comparison  of  Inlet  Velocity  Response  for  Blocked  and  Open 
5  ft,  0.170  in  ID  Lines,  1.0  psi  Terminated  Ramp,  Linear  Results 
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Fig.  68. 
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algorithm  being  used.  The  Initial  flow  field  is  considered  consistent 
with  the  numerical  algorithm  in  this  study  if,  when  introduced  into  the 
algorithm  with  no  signal  present,  continued  iteration  of  the  algorithm 
produces  no  significant  change  in  the  pressure,  velocity,  or 
temperatxire  fields.  The  methoa  of  calculating  this  initial  flow  field 
was  developed  in  Chapter  III  and  reli«ss  on  Eq  (38).  The  initial  flow 
field  is  obtained  with  the  same  friction  model  used  in  the  algorithm 
itself.  This  distinction  is  made  because  empirical  formulas  were  first 
used  to  calculate  the  initial  flow  field  with  very  poor  results. 

Figure  69  shows  the  velocity  profiles  obtained  for  the  initial 
flow  fields  with  the  method  developed  in  Chapter  III.  Figure  69a 
resulted  from  a  pressure  differential  of  0.20  psi  while  Figure  69b  was 
produced  with  a  0.25  psi  differential.  Note  that  the  steady-state 
profile  for  each  flow  regime  has  the  correct  shape.  Figure  69a 
exhibits  the  parabolic  shape  characteristic  of  the  steady  laminar  flow 
velocity  profile.  The  flattened  profile  with  a  steep  gradient  near  the 
wall,  typical  of  turbulent  flow,  is  shown  in  Figure  69b.  These 
qualitative  profiles  are  produced  by  specifying  a  constauit  pressure 
differential  between  the  upstream  and  downstream  ends  of  the 
transmission  line.  Only  the  piecewise  continuous  polynomials  are  used 
to  calculate  the  profile.  The  empirical  and  analytic  expressions 
commonly  used  to  generate  the  steady  velocity  profiles  were  not  used  to 
obtain  these  profiles. 

Tables  8  through  17  are  presented  to  demonstrate  the  performance 
of  the  method  used  to  generate  the  steady-state  profiles  described  in 
Chapter  III.  For  simplicity,  only  six  axial  grid  points 
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and  four  flow  Intervals  were  used  to  produce  these  tables,  however, 
this  method  can  be  used  to  generate  a  steady  flow  field  for  any  number 
of  grid  points  and  flow  intervals. 

Tables  8  through  13  show  the  steady  flow  results  for  air  in  a 
straight,  1.0  ft  long,  0.062  in  ID  tube.  In  Table  8,  the  dimensionless 
pressures  at  each  grid  point  are  shown  with  two  columns  of 
dimensionless  average  velocities.  The  tables  are  in  terms  of  the 
dimensionless  quantities  described  in  Chapter  II  since  this  procedure 
was  designed  to  provide  an  initial  pressure  and  velocity  field  to  the 
numerical  algorithm  in  dimensionless  terms.  The  pressures  in  psia  and 
velocities  in  fps  may  be  recovered  from  the  tables  by  multiplying  the 
dimensionless  pressvures  by  20.538  and  dimensionless  velocities  by 
1128.6,  respectively. 


TABLE  8 

One-Dimensional  Pressures  and  Velocities 


for  Laminar  Flow 

(Ap  =  0.20  psl,  r  =  0.011,  0.022,  0.029,  0.031  in) 


i 

P 

u 

av 

ANALYTIC 

POLYNOKIAL 

0 

0 . 72402376 

0.05523999 

0.05524000 

1 

0.72208675 

0.05538818 

0.05538818 

2 

0.7201444O 

0.05553756 

0.05553756 

3 

0.71P1 

0.0556SG1G 

0.05568817 

4 

0.7x  — .  .1 

0.05584000 

0.05584001 

5 

0.71428571 

0.05599309 

0.05599310 

The  velocities  in  the  column  labeled  "analytic"  were  calculated  using 
the  one-dimensional  compressible  flow  equations  described  in  Chapter 
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IV.  The  values  in  the  column  labeled  "polynomial"  were  calculated 
using  the  piecewise  continuous  polynomial  method  described  in  Chapter 
III.  Values  in  the  column  labeled  "i"  denote  axial  location  along  the 
line  with  the  upstream  end  at  i  «  0  amd  the  downstreaun  end  at 
i  =  5  . 

Wall  shear  stress  values  are  required  to  begin  solving  for  the 
unknown  polynomial  coefficients.  Once  these  coefficients  are  fovind, 
the  velocity  profile  as  well  as  average  velocities  and  velocity 
gradients  cam  then  be  reconstructed.  Close  agreement  between  the 
analytic  and  polynomial  results  for  average  velocities  is  seen  in 
Table  8.  Comparisons  of  the  steady  laminar  velocity  profile  amd 
velocity  gradient  at  each  gridpoint  also  show  close  agreement.  Tables 
9  amd  10  demonstrate  the  accuracy  of  the  piecewise  polynomial 
approximation  for  laminar  flow  at  gridpoint  i  ®  1 

Initial  values  of  the  averaged  velocity  function  in  each  flow 
Interval  are  required  to  begin  the  numerical  solution  of  the  damped 
wave  equations.  The  averaged  velocity  function  produces  the  avei age 
velocity  of  the  flow  in  each  flow  interval.  An  example  of 
dimensionless  values  of  these  functions  for  the  steady  flow  case  are 
shown  in  Table  11.  Moving  from  left  to  right  in  each  row  or, 
physically,  from  the 

center  flow  interval  to  the  wall,  the  decrease  in  the  averaged  velocity 
of  each  flow  interval  can  be  seen. 

Varying  the  locations  of  the  flow  interval  boundaries  has  no 
noticeable  effect  on  the  agreement  between  the  analytic  end 
polynomial  results  for  laminar  flow.  A  different  choice  of  flow 
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TABLE  9 


Comparison  of  Analytic  and  Polynomial  Velocities 


for  Laminar  Flow 


(Ap  =  0.20  psl,  r  =  0.011,  0.022,  0.029,  0.031  In) 

n 


RADIAL 

VELOCITY 

LOCATION 

ANALYTIC 

POLYNOMIAL 

O.OOOOOOOOOOE+00 

1 . 1077636000E-0] 

1. 1077636000E-01 

l.OOOOOOOOOOE-01 

1.0966859640E-01 

1 . 0966859640E-01 

2. OOOOOOOOOOE-01 

1 . 0634530560E-01 

1.0634530560E-01 

3.0000000000E-01 

1 . 0080648760E-01 

1 . 0080648760E-01 

4.0000000000E-01 

9.3052140429E-02 

9.3052142400E-02 

5.0000000000E-01 

8. 3082268 158E-02 

8.30822700C0E-02 

6.0000000000E-01 

7.0896868687E-02 

7.0896870400E-02 

7.0000000000E-01 

5. 64959420 16E-02 

5.6495943600E-02 

8.0000000000E-01 

3.9879491717E-02 

3.9879489600E-02 

9.0000000000E-01 

2. 1047510775E-02 

2. 1047508400E-02 

1 . OOOOOOOOOOE+00 

4.3918815251E-09 

6. 045030750 IE- 13 

TABLE  10 


Comparison  of  Velocity  Gradient 


for  Laminar  Flow 


(Ap  =  0.20  psl,  r  =  0.011,  0.022,  0.029,  0.031  in) 

n 


RADIAL 

LOCATION 

O.OOOOOOOOOOE+00 
l.OOOOOOOOOOE-01 
2.0000000000E-01 
3.0000000000E-01 
4.0000000000E-01 
5.0000000000E-01 
6.0000000000E-01 
7.0000000000E-01 
8.0000000000E-01 
9.0000000000E-01 
1 . OOOOOOOOOOE+00 


SLOPE  OF  RADIAL  VELOCITY  PROFILE 
ANALYTIC  POLYNONIAL 


O.OOOOOOOOOOE+00 
-2.215527Z000E-02 
-4.4310544000E-02 
-6.6465816000E-02 
-8.8621086709E-02 
-1. 1077635871E-01 
-1.3293163071E-01 
-1. 550869027 lE-01 
-1.7724217342E-01 
-1.9939744542E-01 
-2.2155271387E-01 


POLYNONIAL 

0 . OOOOOOOOOOE+00 
-2.2155272000E-02 
-4.4310544000E-02 
-6.6465816000E-02 
-8.S621088000E-02 
-1. 1077636000E-01 
-1.3293163200E-01 
-1.5508690400E-01 
-1.7724217600E-01 
-1.9939744800E-01 
-2.2155272000E-01 


TABLE  11 


Dimensionless  Values  of  Averaged  Velocity  Functions 
for  Laminar  Flow 


(Ap  =  0.20  psi,  r  =  0.011,  0.022,  0.029,  0.031  in) 
n 


FLOW 

INTERVAL 

i 

1 

2 

3 

4 

0 

0.10352470 

0.07570351 

0.0'^431663 

0.00689781 

1 

0. 10380240 

0.07590659 

0.03440868 

0.00691632 

2 

0. 10408236 

0.07611131 

0.03450148 

0.00693497 

3 

0.10436461 

0.07631771 

0.03459504 

0.00695378 

4 

0.10464917 

0.07652579 

0.03468937 

0.00697274 

5 

0.10493608 

0.07673560 

0.03478447 

0.00699185 

interval  boundary  locations  is  shown  in  Tables  12  and  13. 


TABLE  12 

Comparison  of  Radial  Velocity  Profile 
for  Laminar  Flow 


(Ap  =  0.20  psi.  r  =  0.020,  0.024,  0.028.  0.031  in) 

n 


RADIAL 

VELOCITY 

LOCATION 

ANALYTICAL 

POLYNOMIAL 

O.OOOOOOOOOOE+00 

1 . 1077636000E-01 

1 . 1077636000E-01 

1 . OOOOOOOOOOE-01 

1.0966859640E-01 

1 . 0966859640E-01 

2.0000000000E-01 

1.0634530560E-01 

1.0634530560E-01 

3.0000000000E-01 

1.0080648760E-01 

1.0080648760E-01 

4.0000000000E-01 

9.3052142400E-02 

9.3052142400E-02 

5.0000000000E-01 

8.3082270000E-02 

8 . 3082270000E-02 

6.0000000000E-01 

V.0896870400E-02 

7 . 0896870400E-02 

7.0000000000E-01 

b.6495940758E-02 

5 . 6495943600E-02 

8.0000000000E-01 

3.9879492587E-02 

3.9879489600E-02 

9.0000000000E-01 

2. 1047512032E-02 

2. 1047508400E-02 

1 . OOOOOOOOOOE+00 

2.3003039473E-09 

6.0450307501E-13 
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TABLE  13 


Comparison  of  Radial  Velocity  Profile 
for  Laminar  Flow 

(Ap  »  0.20  psi,  r  =  0.020.  0.024,  0.028,  0.031  in) 


SLOPE  OF  RADIAL  VELOCITY  PROFILE 


ANALYTIC 


POLYNOMIAL 


RADIAL 

LOCATION 

0. OOOOOOOOOOE+00 
1 . OOOOOOOOOOE-01 
2.0000000000E-01 
3.0000000000E-01 
4. OOOOOOOOOOE-01 
5.0000000000E-01 
6.0000000000E-01 
7.0000000000E-01 
8.0000000000E-01 
9. OOOOOOOOOOE-01 
1 . OOOOOOOOOOE+00 


0. OOOOOOOOOOE+00 
-2.2155272000E-02 
-4.4310544000E-02 
-6.6465816000E-02 
-8.8621088000E-02 
-1. 1077636000E-01 
-1. 32931 63200E-01 
-1.5508689529E-01 
“1.7724216955E-01 
-1 . 9939744 155E-01 
-2.2155271581E-01 


O.OOOOOOOOOOE+00 
-2.2155272000E-02 
-4.4310544000E-02 
-6.6465816000E-02 
-8.8621088000E-02 
-1. 1077636000E-01 
-1.3293163200E-01 
-1.5508690400E-01 
-1. 77242 17600E-01 
-1 . 9939744800E-01 
-2.2155272000E-01 


Several  other  sets  of  flow  interval  locations  as  well  as  additional 
numbers  of  flow  intervals  were  run  with  the  same  close  agreement 
between  analytical  and  polynomial  results  as  shown  above. 

The  steady  flow  results  shown  above  result  from  a  pressure 
differential  of  0.20  psi  across  the  ends  of  the  tube.  At  this  pressure 
differential,  the  Reynolds  number  based  on  diameter  indicates  laminar 
flow.  When  the  pressure  differential  is  raised  to  0.25,  turbulent  flow 
is  predicted.  Calculation  of  the  velocity  profile  is  no  longer  as 
straightforward  as  it  is  with  laminar  flow.  Agreement  of  the 
analytical  and  polynomial  average  velocities  is  not  nearly  as  good  as 
it  is  in  laminar  flow.  Analytical  and  polynomial  average  velocities 
are  shown  in  Tables  14  through  16.  As  the  pressure  differential  is 
increased,  the  agreement  becomes  better.  It  is  also  possible  to 
improve  the  agreement  at  any  pressure  differential  by  adjusting  the 


TABLE  14 


One-Dlmenslonal  Pressures  amd  Velocities 


for  Turbulent  Flow 

(Ap  =  0.25  psi,  r  »  0.011,  0.022.  0.029,  0.031  in) 

n 


1 

P 

AVERAGE 

VELOCITY 

ANALYTIC 

POLYNOHIAL 

0 

0.72645827 

0.06876869 

0.03337146 

1 

0.72404035 

0.06899834 

0.03349131 

0.72161422 

0.06923032 

0.03360368 

3 

0.71917980 

0.069*6467 

0.03371719 

4 

0.71673699 

0. 069^0142 

0.03383187 

5 

0.71428571 

0.06994062 

0.03395680 

TABLE  15 

One-Dlmensiona] 

Pressures  and  Velocities 

for 

Turbulent  Flow 

3  = 

0.50  psi,  r  = 

n 

0.011,  0.022, 

0.029,  0.031 

1 

p 

AVERAGE 

VELOCITY 

ANALYTIC 

POLYNOMIAL 

0 

0.73863083 

0.07282066 

0.05648626 

1 

0.73382741 

0.07329732 

0.05689787 

2 

0.72899186 

0.07378352 

0.05727394 

3 

0.  •’2412351 

0.07427957 

0.05765761 

4 

0.71922170 

0.07478582 

0.05804913 

5 

0.71428571 

0.07530262 

0.05849654 
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TABLE  16 


One-Dimensional  Pressures  and  Velocities 


for  Turbulent  Flow 

(Ap  =  1.0  psi.  r  *  0.011,  0.022,  0.029,  0.031  in) 

n 


i 

P 

AVERAGE 

VELOCITY 

AMALYTIC 

POLYNOMIAL 

0 

0.76297595 

0.10707231 

0.09229894 

1 

0.75349819 

0.10841910 

0.09365498 

2 

0.74389560 

0. 10981863 

0.09484924 

3 

0.73416304 

0. 11127446 

0.09609093 

4 

0.72429504 

0.11279050 

0.09738332 

5 

0.71428571 

0.11437104 

0.09898636 

locations  of  the  flow  interval  boundaries.  One  possible  reason  for 
disagreement  between  polynomial  and  analytic  results  for  turbulent  flow 
is  the  nature  of  the  turbulent  friction  model  auid  not  the  method  of 
determining  the  steady  velocities.  The  algebraic  turbulence  model  is 
qualitatively  good  since  it  produces  the  correct  profile  shapes. 
However,  a  different  turbulent  friction  model  may  be  required  to 
produce  results  consistent  with  the  average  velocities  produced  by  the 
one-dimensional  empirical  equations  for  steady-flow.  No  work  to 
accomplish  this  has  been  yet  been  done. 

Table  17  is  analogous  to  Table  11  except  the  flow  is  turbulent. 
Comparison  of  the  these  two  tables  shows  the  turbulent  profile 
contained  in  Table  17  is,  as  would  be  expected,  much  flatter  near  the 
center  of  the  line  than  the  laminar  profile  shown  in  Table  11. 


TABLE  17 


Values  of  Averaged  Velocity  Functions 
for  Turbulent  Flow 


(Ap 

*0.25  psi, 

r  =  0.011, 

n 

0.022,  0.029,  0.031  in) 

FLOW  INTERVAL 

1 

1 

2 

3 

4 

0 

0.04493990 

0.04132258 

0.02980270 

0.00327153 

1 

0.04510187 

0.04147166 

0.02990921 

0.00830020 

2 

0.04525317 

0.04161079 

0.03000957 

0.00832807 

3 

0.04540602 

0.04175133 

0.03011096 

0.00835623 

4 

0.04556043 

0.04189331 

0.03021338 

0.00838468 

5 

0.04572928 

0.04204873 

0.03032440 

0.00841455 

The  figures  and  tables  discussed  above  were  all  for  compressible 
flow;  however,  the  routine  used  is  easily  modified  for  incompressible 
flow.  An  initial  flow  field  was  generated  for  a  1.0  in  ID  rigid  line 
filled  with  oil  in  order  to  gauge  the  performance  of  the  numerical 
algorithm  combined  with  the  steady  flow  initialization  process. 
Parameters  used  for  this  system  are  shown  in  Ta*"  18.  Figure  70  shows 
the  transient  response  of  the  system  for  sudden  closure  of  a  valve  at 
the  downstream  end  of  the  rigid  line  with  the  oil  initially  flowing  at 
0.42  fps.  The  initial  velocity  field  was  obtained  using  the  piecewise 
defined  polynomials  described  in  Chapter  III.  At  this  velocity,  the 
Reynolds  number  is  much  less  thain  2300  inuicating  the  initial  flow  is 
laminar.  The  oil  is  treated  as  an  incompressible  fluid  with  a  constant 
sonic  velocity.  The  fluid  velocity  is  less  than  1.0  x  10~*  percent  of 
the  sonic  velocity.  As  a  result,  both  the  linear  and  nonlinear 
numerical  algorithm  produce  the  same  result  since  the  magnitude  of  the 


Table  18 


error  terms  in  the  modified  equation  Eq  (4.39)  is  small.  Data  points 
from  a  paper  by  Zielke  (82:114)  are  shown  in  the  figure  for  comparison. 
Tha  data  points  agree  well  with  the  results  produced  with  either 
algorithm.  Results  from  the  linear  algorithm  are  shown  in  Figure  70. 

It  is  importar  to  note  that  the  finite-difference  model  was  not 
"tuned"  to  fit  the  data.  The  damped-wave  equations  and  numerical 
algorithm  resulted  from  straightforward  manipulation  of  the  basic  fluid 
dynamic  equations  and  polynomial  approximations.  The  use  of  the 
piecewise  continuous  polynomials  in  conjunction  with  the  damped-wave 
equations  appears  to  accurately  represent  this  system. 
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+  -  Zielke,  1968 
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Fig.  70.  Response  of  Oil-Filled  Line  with  0.42  fps  Flow  Velocity  to 
Sudden  Closure  of  Upstream  Valve 


Results  Using  Modal  Superposition 


The  results  presented  for  finite  length  lines  were  obtained  from 
the  finite-difference  solution  of  the  damped  wave  equations  developed 
in  Chapter  II.  The  comparison  of  these  results  with  data  is  good  but 
modest  computer  resources  were  required.  Additi  •)nally,  the  numerical 
algorithm  used  to  solve  the  damped-wave  equations  is  not  compatible 
with  the  time  domain  methods  commonly  used  to  analyze  auid  design 
control  systems.  A  second  type  of  transmission-line  model  which 
provides  a  reasonably  accurate  approximation  to  the  solutions  of  the 
linearized  damped-wave  equations  was  developed  in  Chapter  V. 

The  transient  responses  shown  in  Figure  54  and  in  Figures  58 
through  68  exhibit  periodic  characteristics  related  to  the  physical 
characteristics  of  the  fluid  transmission  line.  The  spikes  in  Figure 
65  and  66  are  separated  by  approximately  2.0  dimensionless  time  units 
throughout  the  response.  The  peaks  of  Figure  54  are  separated  by 
approximately  4.0  dimensionless  time  units  as  the  curve  decays.  This 
periodic  behavior  suggests  natural  frequencies  and  time  constants 
characteristic  of  the  line.  A  model  based  on  the  modal  superposition 
method  should  provide  a  good  basis  for  the  approximation  of  the 
finite-difference  solutions. 

The  transient  response  for  the  blocked,  5  ft,  0.170  in  ID 
transmission  line  filled  with  air  is  shown  in  Figure  71.  Steady-state 
values  for  the  L,  R,  and  C  modal  parameters,  as  developed  by  Lebrun 
(40:459),  wero  used  here.  The  modal  superposition  results  exhibit  the 
periodic  behavior  shown  by  the  linear  finite-difference  results  for  the 
same  transmission  line  configuration  used  in  Figure  54.  Howeve-,  the 
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period  of  the  modal  superposition  results  is  noticeably  different  from 
the  period  of  the  finite-difference  results.  The  waveform  and  the 
damping  also  differ  significantly  between  the  two  models. 

When  unsteady  values  of  L,  R,  and  C  are  calculated  and  used  with 
the  modal  supierposltlon  method  described  in  Chapter  V,  the  result 
greatly  Improves  as  shown  in  Figure  72,  The  diamonds  represent  modal 
supjerpositlon  results  while  the  solid  line  represents  the 
finite-difference  solution  with  laminar  friction  and  no  heat  transfer. 
The  conductance  G  was  set  to  zero  in  Figure  72,  however,  when  values 
for  G  are  calculated  as  described  in  Chapter  V,  the  response  shown  in 
Figure  73  results.  The  modal  superposition  result,  shown  by  the 
triangles,  closely  follows  the  finite-difference  result  although  the 
peaks  are  flatter  and  of  somewhat  smaller  aunplltude  than  in  the 
finite-difference  results.  The  majority  of  the  information  contained 
in  the  waveform  produced  by  the  finite-difference  solution  is  contained 
in  the  waveform  produced  by  the  modal  superposition.  It  is  Impjortant 
to  note  that,  like  the  development  of  the  finite-difference  model, 
modal  superposition  and  all  parameters  used  are  based  on 
straightforward  mathematical  m2Uiipulation  of  the  equations  of  fluid 
dynamics.  No  parameters  have  been  adjusted  to  fit  the  modal  results  to 
the  finite-difference  results.  The  modal  superposition  requires 
significantly  less  computer  resources  than  does  the  solution  of  the 
finite-difference  model.  The  two  modal  superposition  solutions  are 
compared  in  Figure  74.  The  importance  of  the  parameter  G  to  the  modal 
solution  can  clearly  be  seen. 

The  number  of  modes  retained  using  modal  superposition  directly 
affects  the  quality  of  the  solution.  The  minimum  number  of  modes 
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(p-p„)/(p. -p  ) 


TIME  (SECONDS) 


Fig.  71.  1.0  psl  Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Blocked 

Line,  Modal  Superposition  Results 
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Fig.  72.  1.0  psi  Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Blocked 

Line,  Modal  Superposition  Results 
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PRESSURE  (PSIA) 


TIME  (tCo/L) 


Fig.  74.  1.0  psl  Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Blocked 

line,  Compat ison  of  Nodal  Superposition  Results 


required  for  a  desired  solution  accuracy  should  be  retained  to  minimize 
computer  resource  requirements.  The  results  shown  in  Figure  75  were 
obtained  with  only  6  modes.  The  low  quality  of  the  solution  is  evident 
in  the  fluctuations  preceding  the  first  peak  of  the  response.  As  time 
progresses,  the  highest  frequency  modes  are  attenuated  and  become 
increasingly  less  important  in  the  solution.  FigU'‘es  76,  77,  and  78 
involve  10,  15,  and  20  modes  respectively.  The  quality  of  the  solution 
increases  noticeably  with  the  additional  modes.  In  this  case, 
retaining  more  than  20  modes  produced  a  negligible  increase  in  quality. 
The  modal  superposition  result  in  Figure  78  compares  very  well  with 
data  obtained  by  Karam  and  Leonard  (30:501)  for  the  same  system. 

For  the  case  of  the  5  ft  blocked  line  examined  above,  the  modal 
superposition  results  compare  very  well  with  the  finite-difference 
results.  Modal  superpxjsition  must  also  produce  good  results  for  other 
terminations  if  it  is  to  be  as  flexible  as  the  finite-difference  model. 
Figure  79  shows  the  modal  superposition  result  for  a  volume  terminated 
line.  The  modal  sup>erposition  result  tracks  the  finite  difference 
result  quite  closely.  The  modal  superposition  result  contains  the 
inflection  at  the  first  p>eak  exhibited  by  the  finite-difference 
solution.  Although  the  pressure  response  curve  is  fairly  smooth,  the 
velocity  responses  at  the  inlet  of  the  line,  shown  in  Figure  80,  and  at 
the  entrance  into  the  volume,  shown  in  Figure  81,  contain  several  peaks 
and  valleys.  Every  inflection  and  peak  contained  in  the 
finite-difference  solution  is  tracked  by  modal  superposition  in  both 
figures. 

Figure  82  compares  the  response  of  a  branched  system  also  studied 
by  Karam  and  Leonard  (30:501).  The  system  consists  of  a  5  ft,  0.170  in 
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Fig.  75.  1.0  psi  Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Blocked 

Line,  Nodal  Superposition  Resuits  for  6  Modes 
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Fig.  78.  1.0  psi  Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Blocked 

Line,  Modal  Superposition  Results  for  20  Modes 
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Fig.  80.  Transient  /elocity  Response  at  Line  Inlet,  1.0  psi 
Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Volume  Terminated  Line, 
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Fig.  81.  Transient  Velocity  Response  at  Line  Exit  to  Volume,  1.0  psi 


Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Volume  Terminated  Line, 
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ID  and  a  5  ft  0.25  In  ID  line,  both  blocked,  bramching  from  another  5 
ft,  0.170  in  ID  line.  The  response  at  the  blocked  end  of  the  0.25  in 
I.D.  line  is  shown  in  the  figure  with  data  points  from  the  paper  by 
Karaim  and  Leontird.  The  modal  superposition  results  were  calculated 
using  a  separate  modal  superposition  model  and  a  separate  matrix  of  the 
parameters  I,  R,  G,  and  C,  for  each  of  the  three  lines.  Comparison 
with  the  data  points  is  good.  Model  instability  was  experienced  when 
the  inputs  and  outputs  from  each  line  model  were  directly  coupled. 
Instability  occurs  because  the  trunk  line  model  excites  each  branch 
model  at  its  modal  frequencies  since  all  three  lines  are  of  the  same 
length.  This  instability  was  eliminated  by  coupling  the  models  through 
a  volume  equal  to  one-half  of  the  Internal  volume  of  the  5  ft,  0.170  in 
ID  line.  Each  of  the  three  lines  has  approximately  the  same  modal 
frequencies.  Several  values  of  coupling  volume  were  tried  with 
virtually  identical  results.  Coupling  volumes  smaller  than  one  half  of 
the  trunk  line  volume  produced  instability. 

The  results  for  a  single  5  ft  open  line  with  z  1.0  psi  amd  a  0.25 
psi  step  input  are  shown  in  Figures  83  through  86.  As  in  the 
finite-difference  results,  the  solution  "ratchets"  up  to  approach  a 
steady  state.  Figures  83  and  84  show  the  1.0  psi  step  input.  The 
modal  superposition  model  follows  the  finite-difference  model  very 
closely  at  first  but  diverges  significantly  at  later  times.  The 
comparison  between  the  modal  superposition  result  and  the 
finite-difference  result  is  much  better  in  Figures  85  and  86  where  a 
0.25  psi  terminated  ramp  input  signal  is  applied.  The  large 
discrepancy  in  Figures  84  and  85  results  from  the  use  of  the  turbulent 
friction  model  in  the  finite-difference  model.  The  values  of  I,  R,  G, 
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Fig.  82.  Transient  Response  in  One  Leg  of  Branched  Line,  1.0  psi 
Terminated  Ramp  Input  in  5  ft,  0.170  in  ID,  Blocked  Branches,  Line  2 


185 


and  C  used  in  the  modal  superposition  model  were  calculated  using 
laminar  friction  only.  No  method  of  determining  values  for  these 
parameters  for  unsteady  turbulent  friction  has  yet  been  developed  for 
this  application.  As  stated  earlier,  methods  have  been  developed  to 
incorporate  turbulent  friction  resulting  from  an  established  mean 
turbulent  flow  into  these  parameters  when  determining  frequency 
response.  Modal  superp>oslticn,  however,  requires  a  method  to  account 
for  the  unsteady  effects  of  transition  to  turbulence  in  an  initially 
undisturbed  flow  for  a  signal  of  arbitrary  shape.  This  poses  a 
difficult  problem. 

The  results  from  modal  superposition  demonstrate  a  capability  to 
accurately  predict  the  linear  transient  response  for  laminar  flow. 

Much  more  work  is  required  to  give  the  modal  superposition  model  the 
ability  to  track  the  nonlinear  finite-difference  results,  as  well  as 
account  for  effects  due  to  unsteady  turbulent  flow.  Modal 
superposition,  although  only  strictly  valid  for  linear  systems,  has 
been  modified  to  model  nonlinear  transient  responses  in  support  cable 
systems  (47).  Methods  used  in  these  cable  models  may  be  applicable  to 
transmission  lines. 
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VELOCITY  AT  X  =  0  (EPS) 


Fig.  83.  Transient  Velocity  Response  at  Line  Inlet,  1.0  psi 
Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Open  Line 
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Fig.  84.  Transient  Velocity  Response  at  Line  Exit,  1.0  psl 


Temlnated  Ramp  Input  in  5  ft,  0.170  in  ID  Open  Line 
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Fig.  85.  Transient  Velocity  Response  at  Line  Inlet,  0.25  psi 
Terminated  Ramp  Input  in  5  ft,  0.170  in  ID  Open  Line 
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Fig.  86.  Transient  Velocity  Response  at  Line  Exit,  0.25  psi 
Teminated  Ramp  Input  in  5  ft,  0.170  in  ID  Open  Line 


VII.  Conclusions  and  Suggestions  for  Further  Research 


Conclusions 


A  finite-difference  model  and  a  modal  superposition  model  capable 
of  predicting  the  waveforms  of  flow  and  pressure  transients  from  input 
signals  of  arbitrary  shape  aind  finite  amplitude  were  presented  in  this 
study.  The  finite-difference  model  improved  on  the  approach  of  Tsao 
(67)  by  using  the  numerical  algorithm  developed  by  Gabutti  (21)  to  more 
correctly  represent  the  physics  of  the  problem.  A  method  to 
algebraically  generate  a  steady-state  velocity  field  was  also  developed 
to  provide  initial  conditions  for  cases  where  a  steady,  mean  flow  is 
present.  A  simple,  turbulent-friction  model,  modified  for  transient 
flow,  was  included. 

The  modal  superposition  model,  based  on  the  work  of  Lebrun  (39), 
was  shown  to  closely  track  the  important  features  of  the  linear, 
finite-difference  results  when  modified  to  use  parameters  based  on 
well-known,  frequency-domain,  transmission-line  theory. 

The  finite-difference  model  provides  the  most  general  method  by 
which  to  predict  the  behavior  of  transient  flow  in  a  fluid  transmission 
line.  The  finite-difference  model  incorporates  unsteady  friction  amd 
heat  transfer,  both  laminar  and  turbulent,  as  well  as  the  effects  of  an 
established  mean  throughflow.  This  approach  is  similar  to  that 
presented  by  Tsao  (69),  in  that  it  uses  piecewise  defined  polynomials 
to  reduce  the  number  of  grid  points  needed  in  the  radial  direction. 

The  transient  velocity  profile  at  each  axial  station  on  the 
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computational  grid  can  be  predicted  with  relatively  modest  computer 
resources. 

Two  improvements  to  Tsao’ s  approach  have  been  made  in  this  study. 
An  algorithm  developed  by  Gabutti  (21)  to  represent  correctly  the 
physics  of  a  nonlinear  problem  was  applied  to  a  set  of  nonlinear, 
damp>ed-wave  eqxiatlons.  For  the  linear  case  alone,  many  numerical 
schemes,  such  as  the  Lax-Wendroff  algorithm,  will  correctly  represent 
the  physics  of  the  problem  when  used  to  solve  the  damped 
wave-equations.  When  nonlinear  effects  are  to  be  considered,  more 
attention  to  the  physics  of  the  problem  is  required  to  select  the 
proper  algorithm.  Many  grid  points  are  required  to  model  the  progress 
of  a  signal  traveling  through  a  transmission  line.  Failure  to  account 
correctly  for  the  physics  of  the  problem  produces  errors  at  each  grid 
point  which  are  propagated  along  with  the  signal  and  are  cumulative. 

Use  of  the  Gabutti  algorithm  reduces  this  source  of  error  as  well  as 
providing  a  logical  aind  physically  consistent  meauis  to  handle  boundary 
conditions  present  at  the  ends  of  the  transmission  line. 

The  second  improvement  is  the  use  of  the  piecewise  defined 
polynomials  to  generate  a  steady-state  velocity  field  to  model  a 
transmission  line  with  a  steady,  laminar  or  turbulent,  meam  flow.  A 
relationship  between  the  we 11 -developed,  one-dimensional,  steady-flow 
equations  amd  piecewise  defined  polynomials  was  developed  for 
steady-flow  velocity  profiles  recognizing  that  no  net  momentum  transfer 
occurs  across  the  flow  interval  boundaries  during  steady  flow.  The 
initial  flow  field  for  a  transmission  line  with  a  steady,  mean  flow  is 
now  obtained  without  carrying  the  solution  of  the  damped-wave  equations 
out  to  steady  state. 
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The  finite-difference  model  allows  the  use  of  separate  models  to 
predict  the  effects  of  unsteady  friction  and  heat  transfer.  Studies  by 
Manning  and  Sandoz  used  an  experimentally  determined  friction  factor  to 
match  a  one-dimensional  model  to  experimental  results.  No  separate 
provision  was  made  for  heat  transfer,  the  effect  of  which  was  assumed 
to  be  lumped  into  the  friction  factor.  The  finite-difference  model  in 
this  study  was  shown  to  produce  results  consistent  with  published 
experimental  results  for  finite-length  lines  without  the  need  for  trial 
and  error  determination  of  damping  coefficients.  The  finite-difference 
model  developed  in  this  study  is  superior  to  the  Sandoz  and  Manning 
(60)  model  in  this  respect. 

The  finite-difference  model  is  limited  by  the  amplitude  of  the 
input  signal.  The  model  was  developed  from  the  nonconservative  form  of 
the  Navier-Stokes  equations.  As  a  result,  this  model  fails  in  the 
presence  of  discontinuities.  Just  like  the  method  of  characteristics. 
The  performance  of  this  model  depends  on  amplitude  of  the  input  signal, 
the  input  waveform,  transmission- line  diameter,  and  the  length  of  line 
over  which  signal  must  travel.  This  model  was  demonstrated  in  a 
0.305  in  ID  semi-infinite  line.  For  the  positive  half-sinusoidal 
signal  of  10  psi  amplitude,  a  discontinuity  formed  in  the  solution  at  a 
distance  of  10  ft.  Discontinuities  form  at  greater  distances  for 
signals  of  smaller  amplitudes.  Sufficiently  small  signals  will  never 
produce  solutions  with  discontinuities  since  the  frictional  and  heat 
transfer  damping  terms  contained  in  the  model  attenuate  the  signal 
below  amplitudes  required  to  form  discontinuities.  Signals  with  steep 
slopes,  such  as  steps  or  square  pulses,  will  produce  discontinuities 
much  faster  than  half-sinusoidal  signals.  Also,  lines  of  smaller 
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diameter  produce  greater  frictional  damping,  delaying  formation  of 
discontinuities  in  the  solution  by  attenuating  signals  more  rapidly 
thain  larger  lines. 

Results  from  the  finite-difference  model  demonstrate  that  the 
threshold  for  nonlinear  behavior  is  much  lower  than  1.0  psi.  This  is 
in  agreement  with  Sandoz  and  Manning  who  say  that  1.0  psi  is  well 
within  range  of -nonlinearity  (58).  Significant  nonlinear  behavior  was 
shown  at  input  signal  amplitudes  of  0.25  psi.  No  signals  with 
amplitudes  below  0.25  psi  were  tried.  However,  linear  and  nonlinear 
results  converge  as  the  aimplitude  of  the  Input  signal  decreases. 

The  unsteady  turbulence  model  used  in  the  finite-difference  method 
provides  a  qualitative  idea  of  the  effects  of  unsteady  turbulent 
friction,  but  is  very  rudimentary.  Much  more  data  is  needed  to  improve 
the  turbulent  friction  model.  Improvements  are  needed  to  produce 
consistent  results  for  steady  flows  with  auiy  location  of  flow  Interval 
boundaries.  Poor  performance  in  predicting  steady-flow  velocity 
profiles  does  not  invalidate  the  results  for  transient,  turbulent  flow. 
Most  of  the  damping  in  transient  flow  occurs  near  the  wall,  away  from 
the  outer  layer  where  the  eddy  viscosity  model  may  not  be  valid.  Also, 
much  more  data  is  needed  on  the  transition  of  unsteady  laminar  flow  to 
unsteady  turbulent  flow.  No  clues  as  to  what  significant  parameters  to 
base  transition  on  were  found  in  the  literature  during  this  study. 

Some  previous  studies  have  experimentally  observed  that  damping 
significantly  increases,  with  signal  amplitude  in  long,  coiled  lines 
(27:280;  29:500-501).  Braden  (5)  developed  a  model  that  exhibited  this 
increase  in  damping  with  signal  amplitude.  However,  Braden 
specifically  designed  his  model  to  produce  this  effect.  No  results 
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matching  those  of  the  previous  studies  were  observed.  No  mechanism  was 
included  in  the  finite-difference  model  to  mimic  these  results  because 
the  source  of  this  aunpHtude  dependent  damping  observed  in  the  previous 
studies  is  not  well  understood.  The  nonliner.r  damping  effect  observed 
may  result  from  the  effects  of  the  colled  lines  or  an  effect  of 
transition  to  turbulence  in  the  transient  flow,  from  an  initially 
undisturbed  state.  As  can  be  seen,  development  of  unsteady  transition 
parameters  for  use  with  this  model  promises  to  be  a  difficult  problem. 
However,  in  spite  of  its  weaknesses,  this  finite-difference  model  is  a 
suitable  tool  for  designing  experiments  to  better  characterize  fast 
transient  flow. 

The  second  model  developed  in  this  study  was  the  modal 
superposition  model.  This  model  produced  results  consistent  with 
results  from  the  f inlte-dlfference  model.  Modal  superposition  as 
developed  in  this  study,  however,  is  limited  to  linear  systems  with  no 
established  mean  flow.  This  model  was  included  in  this  study  as  a 
means  to  approximate  the  results  of  the  finite-difference  model 
discussed  above.  The  modal  superposition  model  is  convenient  to  use  in 
the  design  of  fluid  networks  and  requires  much  less  computational  power 
than  the  finite-difference  model.  Modifications  to  modal  superposition 
models  allowing  the  approximation  of  nonlinear  behavior  have  been 
successful  in  structural  mechanics.  Similar  modifications  to  the  model 
developed  in  this  study  may  be  possible. 

The  strategy  is  to  use  the  modal  superposition  model  and 
finite-difference  model  as  a  pair.  The  finite-difference  model,  itself 
too  cumbersome  for  fluid  system  design,  is  used  to  characterize 
individual  transmission  lines  and  signals  being  considered  for  use  in  a 
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fluid  system.  The  modal  superposition  model  is  then  adjusted  to  model 
finite-difference  results  for  individual  transmission  lines  and,  in 
turn,  is  used  to  model  entire  fluid  systems. 

Modal  superposition,  as  shown  in  the  previous  section  of  this 
study,  closely  tracked  the  important  features  of  the  linear 
finite-difference  results  for  finite  lines.  The  model  and  all 
parameters  used  were  obtained  by  straightforward  manipulatior  of  the 
governing  fluid  dynamic  equations  with  no  trial  and  error  Involved. 
Model  parameters  based  on  the  well  known  linear  transmission  line 
theory  were  used  to  obtain  transient  responses  directly  without  the 
need  to  invert  from  the  frequency  domain  to  the  time  domain.  Proper 
modification  of  the  I,  R,  C,  and  G  parameters  used  in  this  model  would 
allow  modal  superposition  to  simulate  some  of  the  nonlinear  effects 
observed  in  the  finite-difference  results.  This  dual  model  approach 
may  provide  a  useful  design  tool  when  nonlinear  effects  must  be 
considered  in  fluid  syste..*.  design. 

Suggestions  for  Further  Research 

The  results  for  the  finite-difference  model  presented  in  Chapter 
VI  show  that  the  frictional  and  heat  transfer  damping  functions  are 
critical  to  accurate  prediction  of  the  transient  response  of  a  fluid 
transmission  line  to  an  applied  signal.  In  the  ca^e  of  laminar 
friction  and  heat  transfer,  evaluation  of  these  damping  functions  is 
straightforward,  however,  transient  turbulent  friction  and  heat 
transfer  are  poorly  understood.  With  the  exception  of  shock  tube 
studies,  no  information  regarding  the  transition  to  turbulence  of  an 
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initially  undisturbed  fluid  was  found  during  the  course  of  this  study. 
Experimental  work  examining  this  type  of  transition  for  single  pulse 
and  step  input  signals  must  be  done  before  an  accurate  representation 
of  this  effect  can  be  incorporated  into  a  finite-difference  model. 

Data  is  needed  to  understand  how  to  properly  model  unsteady  turbulent 
friction.  A  more  sophisticated  model  such  as  the  model  of  Nee  and 
Kovasznay  may  be  required  (46).  This  model  consists  of  a  single 
partial-differential  equation  in  which  eddy  viscosity  is  treated  as  a 
conserved  property.  Other,  more  recent  models  may  also  be  suitable. 

This  experimental  work  may  not  have  been  done  due  to  the  extreme 
difficulty  in  generating  pulses  of  sufficient  amplitude  amd  the 
desired  shape.  Some  experimental  work  was  originally  intended  as  part 
of  this  study  but  the  difficulties  that  arose  in  creating  the  required 
input  signals  prevented  the  production  of  any  useful  results.  Figure 
87  shows  a  diagram  of  the  experimental  setup  that  was  tried. 


Fig.  87.  Transmission  Line  Transient  Measurement  Instrumentation 
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The  signal  generator  consisted  of  two  quick  acting  electrically 
operated  solenoid  valves.  Signals  produced  by  this  means  proved  to  be 
of  Insufficient  amplitude  and  poor  shape.  A  number  of  generators  for 
sinusoidal  signals  are  available  for  small-amplitude  signals.  The 
development  of  generators  capable  of  reliably  generating  well-shaped, 
high-aunplitude  pulses  is  a  requirement  for  continued  work  in  transient 
turbulent  friction  in  transmission  lines. 

Much  additional  work  is  necessary  to  incorporate  the  turbulent 
friction  and  heat  transfer  effects  as  well  as  the  nonlinear  effects 
produced  by  the  finite-amplitude  signals  into  the  modal  model 
introduced  in  this  study.  As  described  in  Chapter  V,  the  modal 
parameters  for  this  model  are  calculated  from  approximations  of 
frequency-domain  expressions  for  laminar  flow.  Modifications  to  these 
frequency-domain  expressions  for  I,  R,  C,  and  G  are  needed  to  account 
for  the  effects  of  turbulence.  Some  modifications  to  these  expressions 
are  available  for  turbulence  resulting  from  the  presence  of  am 
established  meam  turbulent  flow  but  nothing  presently  exists  to  include 
the  effect  of  transient  turbulence.  To  the  best  of  this  author’s 
knowledge,  this  study  is  the  first  direct  use  of  the  frequency  domain 
expressions  as  parameters  in  a  model  of  this  type.  The  results  of  this 
model  are  good  enough  to  Justify  further  development  of  this  approach. 

Incorporation  of  the  nonlinear  waveform  evolution  indicated  by  the 
finite-difference  model  into  the  modal  superposition  model  also 
remains.  In  addition,  a  means  to  allow  this  model  to  interact  with  a 
meam  flow  is  necessary.  As  it  exists  now,  modal  superposition  is 
limited  to  small  amplitude  signals  traveling  in  an  undisturbed  fluid. 
Two  techniques  which  might  be  applied  to  the  modal  model  presented  in 


this  study  have  been  published  by  Moriss  (45)  and  Thompson  and  Kohr 

(66). 

Moriss  (45)  applied  the  modal  superposition  method  to  the 
calculation  of  the  nonlinear  dynamic  response  of  structures.  The 
dynamic  equations  discussed  by  Moriss  are  similar  in  form  to  those  of 
the  modal  model  of  this  study.  The  nonlinear  effects  are  introduced  by 
modifying  the  stiffness  matrix  with  functions  that  change  with  the 
deformed  position  of  each  structural  member  (45  :66).  Use  of  this 
technique  for  the  transmission  line  modal  model  would  require  functions 
to  modify  expressions  for  I,  R,  C,  and  G,  presently  calculated  at  the 
beginning  of  the  solution.  Frequency-dependent  expressions  for  these 
parameters,  valid  for  finite-amplitude  signals  in  fluid  transmission 
lines,  have  been  presented,  however,  their  incorporation  into  modal 
superposition  is  not  straightforward.  Functions  which  change  with 
pressure  might  be  one  alternative.  Another  potential  technique  uses 
sensitivity  analysis  to  develop  sensitivity  functions  which  are  then 
minimized  with  respect  to  the  desired  system  performamce  curve  (66). 
This  technique  was  developed  for  use  in  compensating  nonlinear  control 
systems.  This  approach  offers  a  means  to  tune  the  modal  superposition 
model  to  mimic  results  from  the  finite-difference  model. 

The  modifications  discussed  above  do  not  address  interaction  with 
a  mean  flow.  Methods  used  to  model  acoustic  streaming  and  propagation 
of  periodic  signals  of  finite  amplitude  may  be  of  use  here.  Much  work 
has  been  done  in  the  area  of  nonlinear  acoustic  waveguides  for 
sinusoidal  signals.  A  number  of  papers  in  this  area  were  encountered 
in  the  search  for  information  relating  to  the  propagation  of 
finite-amplitude  signals.  Further  examination  of  literature  in  the 
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area  of  nonlinear  acoustics  should  provide  a  basis  for  modification  of 
modal  superposition  to  interact  with  mean  flow  and,  possibly,  a  means 
to  modify  the  stiffness  matrix  of  the  modal  model  to  account  for 
nonlinear  behavior. 

The  task  of  further  development  of  this  model  will  not  be  simple 
and  much  more  work  is  necessary  before  this  model  does  everythirig 
discussed  above. 
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APPENDIX  A 


Transient  Solutions  Using  The  Brown-Nlchols  Model 


An  analogy  drawn  between  linear  electrical  and  fluid  transmission 
lines  of  uniform  diameter  Is  useful  In  obtaining  the  transient  response 
of  a  flu^d  traiisffiijslcn  I’.ne.  The  linear  dynamic  characteristics  of  a 
fluid  transmission  line  may  described  by  four  distributed  parameters 
(.7:20).  These  parameters  are  I  (Inductamce,  also  known  as  Inertance), 
R  (resistance),  G  (shunt  conductance),  and  C  (shunt  capacity  or 
capacitance).  These  four  parameters  are  combined  to  produce,  Z  =  R+JI 
(Impedance)  amd  Y  ®  G+JC  (admlttauice) .  Z  and  Y  are  also  combined  to 
produce  F  (propagation  operator)  and  Z  (characteristic  Impedance). 

C 

The  propagation  operator  and  characteristic  Impedance  result  from  the 
solution  of  the  differential  equations  that  arise  when  R,  I,  G,  and  C 
are  assumed  to  be  distributed  along  the  length  of  a  uniform  fluid 
transmission  line  as  shown  In  Figure  A. 1.  The  transmission  matrix  Is 
obtained  from  Figure  A. 1  in  terms  of  the  Laplace  variable  s  (36:67): 


P 

1 


cosh  r(s)x 
2  sinh  r(s)x 


Z  sinh  r(s)x 

C 

cosh  r(s)x 


(A.  1) 


where  the  subscripts  1  and  2  denote  values  at  0  and  x,  respectively, 
and  P  and  Q  are  pressure  and  flow,  respectively.  r(s)  and  Z  determine 

C 

the  response  of  the  transmission  line. 
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Fig.  A.l.  Equivalent  Circuit  of  an  Infinitesimal  Portion  of  Uniform 
Trauismission  Line  in  Time  Domain 


The  momentum,  continuity,  energy,  and  state  equations,  simplified 
for  a  fluid  transmission  line  of  uniform  diameter,  are  written  (22:4) 
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■  a^u  1  au  ■ 
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ap  r  a^T  1  aT 

—  =  a  — ^  + - 

at  ar  r  ar 


dp  dp 

—  =  —  (for  liquids) 


dp  dT  dp 

—  +  —  =  —  (for  gases) 


(A. 2) 


(A. 3) 


(A. 4) 


(A. 5a) 


(A. 5b) 
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These  equations  are  solved  in  the  Laplace  domain  in  conjunction  with 
the  expressions 


Z(s)  F(s)  =  - 


Y(s)  P(s)  =  - 


ap(s) 

ds 

aF(s) 

as 


F(s)  =  2np  I  rV(s)dr 
ojo 


to  produce  expressions  for  Z(s)  and  Y(s)  (7:549) 


Z(s)  = 


I  s 

a 


1  + 


2J^ (JaVs/v^) 


AaVs/v  J  {.AaVs/v  ) 
0  0  0 


Y(s)  =  C  s 


1  + 


ZCy-l ) ( Ja/Prs/v^ ) 


JavPrsTv  J  (JaVPrsTv  ) 
0  0  0 


P  A 

where  I  =  -  auid  C  =  -  Equations  (A.  7)  and  (A.  8)  can  be 


a  2 

pa 


written  in  the  frequency  domain  to  produce  (49:7) 


Z(w)  =  jwl 


1  - 


Y(u)  =  JoC  \ 


1  + 


2Jj 

2(y-l)Jj  ^J^^V8Prw/w^  j 


(A. 6a) 

(A. 6b) 

(A. 6c) 

(A. 7) 

(A. 8) 


(A. 9) 

(A. 10) 
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The  parameter  pairs: 


u 

Sirv 

=  —7—;  (i)  = 

8irv 

(A. 11a) 

V 

A  T 

A 

a 

V 

=  V2v/(i);  5^ 

=  VZv^/u) 

(A. 11b) 

D/2 

D/2 

h 

S  ^  S 

(A. 11c) 

a  ’  " 

V 

a 

T 

may  be  used  to  simplify  Eqs  (A. 9)  and  (A. 10)  (49:6-7): 


ZCu)  =  1 


r 


Y(u)  =  JwCj  1 


V 


2J,  ^  hj 

^  h^' 


(A. 12) 


(A. 13) 


For  a  semi-infinite  line,  Brown  inverts  the  expression  p  =  e 

to  obtain  the  transient  response  for  a  semi-infinite  line  (7:551).  For 
a  terminated  transmission  line  of  length  L,  a  transfer  function 
relating  to  is  obtained  using  Eq  (A.l): 


P  1 

—  =  = -  =  G  (s)  (A. 14) 

P 

1  —  sinh  FL  +  cosh  FL 

Z 

L 


This  transfer  function  forms  the  basis  of  the  frequency  response 


conversion  and  the  fast  Fourier  transform  methods  for  obtaining  the 
transient  response  of  a  terminated  transmission  line. 

The  frequency  response  conversion  method  uses  Eq  (A. 14)  with  s 
replaced  by  Jw.  The  input  signal,  P^(t)  is  represented  by  a  Fourier 
series,  P  («  t).  G(Jw)  is  written  in  polar  form  as  G(w  )  -  M  for 

Ik  kk  k 

each  discrete  frequency  resulting  from  the  Fourier  series.  The 
expression  for  the  transient  response  is  written  (33:4) 

« 

P^(t)  =  y  P, t  +  d  )  (A. 15) 

2  L  k  1  k  k 

ksO 

A  related  method  is  the  fast  Fourier  transform  method.  The  discrete 
Fourier  transform  of  the  input  signal  is  computed  by 

M-l 

Pi(j  =  At  J]  P^(kAt)e’-’^”'''‘^'*  (A.  16) 

k=0 

For  ,  the  Fourier  transform  of  the  output  is  described  by 

P^Ok)  =  GOj^)  Pj(JUj^)  (A. 17) 

The  output  waveform  P  (t)  is  obtained  from  the  inverse  transform  of  the 

2 

complex  conjugate  of  Eq  (A. 17)  (71:214): 

N-l 

P^OcAt)  =  Af  J]  P^(J27mAf)e"-’^"'‘''"  (A.  18) 

n  =  0 
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A  third  method  for  obtaining  the  transient  solution  relies  on  a 
modal  approximation  of  the  terms  of  Eq  (A.l)  leading  to  a  state  space 
representation  (24;  25).  The  state  variables  are  chosen  to  be  (P  ) 

2  n 

and  (Z  Q  )  where  n=l,2,  •  •  ‘N  and: 
c  1  n 


Z  Q  -  (Z  Q  ),  *  (Z  Q,  ),  + 

cl  cll  cl2 


(A. 19a) 
(A. 19b) 


The  dynamics  of  the  transmission  line  are  described  by  the  expression 


z  3  A  X  -«■  B  u 


(A. 20) 


where  the  components  of  matrices  A  and  B  are,  themselves,  matrices: 


A  = 

r  A  1 

'  A 

2 

where  A  = 

fA 

11 

A 

12 

i 

A 

A 

A 

21 

22 

(A. 21a) 


2  where  B  = 


B  B 

11  12 

B  B 

21  22 


(A. 21b) 


(P) 

2  n 
0  1  o 


(A. 21c) 
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Calculation  of  the  scalar  values  contained  in  the  A  and  B  matrices  is 
algebraically  complex  (25). 

An  approach  based  on  the  method  of  characteristics  auid  Eq  (A. 7) 
was  also  developed  (9:219): 


1  dp  dV  f  Z 

o 

± - +  —  +  —  ±aa 

al  dt  dt  [  I 


(A. 22) 


This  approach  is  a  generalization  of  the  one  presented  by  Zlelke  (82). 
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APPENDIX  B 


Derivation  of  Damped-Wave  Equations 


The  continuity,  momentum,  energy,  auid  state  equations  in  terms  of 
cylindrical  coordinates  with  mass  averaged  dependent  variables  are 
written  as  follows  (58:11-12): 

Continuity 


a(rp)  a(rpu)  3(rpv) 

-  +  -  +  -  =  0  (B.l) 

at  ax  dr 


x-momentum 


a(rpu) 


at 


puv)' 


=  0 


(B.2) 


r-momentxam 


a 

—  (rpv) 

at 


//  // 
puv 


+ 


+  pvv 


+  T 


96 


P 


//// 

-  pww  = 


0 


(B.3) 
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Total  Energy 


a  ^  a 

— (re)  +  — r 

at  ax 


-  -  -  - 

u(e  +  p)  +  q  +  puh  -  u(t  -  puu)  - 

X  XX 


v(t  -  pvu) 

rx 

M 

-  u 

- im 

puu 

^xx  ■  2 

-  V 

Mil 

(  \ 

pvu 

^rx  ‘  2‘ 

k  j 

+  — r 
ar 


v(e  +  p)  +  q  +  pvh  - 


_  it  tk  ^  titt  $i 

u(t  -  puv)  -  v(t  -  pw)  -  U 
xr  rr 


II II  ^ 

puv 


-im 


{  _  pw' 

M 

-  W 

(  \  1 
pwv 

rr  2 

0r  2 

k  ^ 

)  - 

=  0  (B.4) 


Equation  of  State 


dp  dp  dT 

p  p  T 


(B.5) 


The  viscous  layer  formed  by  transient  flow  in  a  long  fluid 
line  is  thin  with  respect  to  the  line  length  provided  that 
^  «  1  If  the  largest  viscous  term  is  assumed  to  be  of  approximately 
the  same  magnitude  as  the  inertia  terms,  Eqs  (B.l)  through  (B.5)  can 
be  greatly  simplified  by  applying  the  boundary  layer  assumptions 
(2:208).  These  boundary  layer  assumptions  allow  the  second  derivatives 
in  the  direction  of  the  flow  to  be  neglected  compared  to  corresponding 
derivatives  transverse  to  the  flow.  In  other  words,  the  heat  transfer 
and  momentum  transfer  in  the  radial  direction  are  significantly  greater 
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than  in  the  axial  direction  for  a  long,  cylindrical  line.  The  boundary 
layer  assumptions  also  allow  the  momentum  equation  for  flow  in  the 
radial  direction  to  be  neglected. 

This  simplification  process  is  done  methodically  with  an  order  of 
magnitude  analysis  by  assuming  gradients  perpendicular  to  the  surface 
are  of  0(1)  and  gradients  parallel  to  surface  are  0(5)  (14:66).  Also, 

_  ^  It  H  HM 

for  1  V  Mach  number;  t  »  t  and  puv  »  puu  (14:74).  Although 

xr  XX 

pulses  travel  at  the  local  speed  of  sound,  the  local  fluid  velocities 
induced  by  the  passage  of  these  signals  are  much  less  than  Mach  1 .  As 
a  result,  only  the  x-momentum  equation  for  the  axial  velocity  component 
u  is  required.  The  variable  u  is  a  function  of  t,  x,  and  r  where  r  is 
measured  from  the  centerline. 

Despite  the  assumption  of  parallel  flow,  the  turbulence  remains 

ft  / 

three  dimensional  with  fluctuating  velocity  components  u  ,  v  ,  and  w 
in  the  x,  r,  and  9  directions,  respectively.  Based  on  the  above 
assumptions,  the  x-momentum  equation  Eq  (B. 1)  may  be  written  as 
follows: 


3  a  ___  _ 

— (rpu)  +  — [r(puu  +  p)] 

at  dx 


+ 


ft  If 

puv)] 


=  0 


(B.6) 


These  same  assumptions  allow  the  total  thermal  energy  equation,  Eq 
(B.4),  to  be  written  as 


a  (re)  a 

-  +  —  r[u(e  +  p)  -  V 

at  ax 

a 

—  r[v(e 
ar 


_  //  ti 

(t  -  pvu)l  + 

rx  J 


—  —  " 

-*•  p)  +  q  +  pvh  -  u(t 

r  xr 


//  n 

puv ) ]  =  0 


(B.7) 
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Terms  such  as  u(t  -  in  Eq  (B.4)  are  usually  neglected  since  they 

xr  2  _ 

^  ^  H  II 

are  much  smaller  them  thermal  dissipation  terms  like  u(t  -  puv) 

xr 

(57:5). 

Under  the  boundary  layer  approximation,  u  =  u  but  v  *  v 
(14:75).  Even  though  parallel  flow  was  assumed,  v  is  not  zero.  The 
definition  of  the  mass  weighted  variable  allows  v  to  be  written  as 


pv  (p  V  +  pv) 


(B.8) 


With  V  equal  to  zero  from  the  boundary  layer  assumptions,  Eq  (B.8) 


reduces  to  v  =  —  .  In  a  thin  shear  layer,  pv  is  of  0(5)  as  compared 
P 

_  ^  _  “7777 

to  u  which  is  of  0(1)  (14:75).  This  allows  the  terms  v  (t  -  pvu)  and 

PX 

v(e  +  p)  to  be  neglected  in  relation  to  the  terms  u(e  +  p)  and 

^  ////  ~  H  tt 

q^  +  pvh  -  u(t^^-  puv)  in  the  energy  equation.  The  terms  retained  in 
the  energy  equation  are  all  of  order  0(1).  In  the  momentum  equation, 

—  ////  ^  ft  n 

the  term  pvu  -  t  +  puv  reduces  to  -  t  +  puv  .  It  is  now  possible 

xr  xr 

to  write  Eqs  (B.6)  and  (B.7)  solely  in  terms  of  time  averaged 
variables.  The  momentum  equation  may  be  written  as 


a  __  a  _  _  d  _-r, 

—  (rpu)  +  — [r(puu  +  p)]  +  — [r(puv  -  t  )]  =  0  (B.9) 


Note  that  the  term  puv  in  Eq  (B.6)  has  been  changed  to  puv  in  Eq  (B.9) 
(14:73). 


Since  e  =  e  under  the  boundary  layer  approximations, 
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_  _  -2 

e=ph-p+p^  (14:75).  Expansion  of  e  in  Eq  (B.7)  gives  the 
following  version  of  the  energy  equation: 


3t 


ph-p 


— 2^  a 

pu 


0  r  f 

—  ~  Q\X 

+  — r  u  ph  +  + 

ax 

L  V  /  J 


A  mm  *  ' 

— r[  q  +  pvh  -  u  (t  -  puv)]  =  0 


(B.IO) 


ar 


The  energy  equation  can  be  put  into  a  more  convenient  form  by 
subtracting  the  mean  kinetic  energy  equation  from  the  total  energy 
equation  where  the  mean  kinetic  energy  equation  is  written  as  follows: 


a 

a 

f  ““A 

*  O 

+  — 

-pu 

ru% 

at 

ax 

ap  a 


+  ru - r[u(T  -  puv)]  +  r(T  -  puv)—  =  0 

ax  ar  '■*  ar 


(B.ll) 


The  resulting  energy  equation  is  written  below: 


^  __  a(rp)  a  _  a  _  _-r,  _ap  _  -~'du 

—  (rph)  -  -  +  — (ruph)  +  — r(q  +  pvh)  +  ru —  -  r(T  -  puv) —  =  0 

at  at  ax  ar  ''  ax  ar 

(B.12) 

_—r,  au 

The  term  (t  -  puv) —  is  the  dissipation.  Cebeci  auid  Smith  show 

r  X 

ar 

a  nondlmensional  form  of  this  equation  in  nonconservative  form  as 


,ah  ,  *ah 
p  —  +  pu  —  =  Ec 

at  ax 


'  •  •' 

ap  ,ap 

1 

a 

p*ah*  .  -r,* 

—  +  u  - 

+  — 

— 

r - -  Rl  p  (hv) 

at  ax 

\  ) 

Rl 

ar 

_  Pr  ar 

Ec 

Rl 


.5u  ,  — / 

p  —  -  Rl  p  (uv) 

ar 


au 

ar 


(B. 13) 
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The  asterisk  denotes  nond 1 mens Iona 1  quantities  (14:77).  The  Eckert 

2 

U 

number  for  pipe  flow  is  defined  as  E  >  h  while  the  Reynolds 

h 

UcPL  “ 

number  is  defined  as  R  =  — —  .  u  denotes  centerline  velocity.  Low 

L  c 

fluid  velocities  produce  a  small  E  while  R  is  large  within  the 
boundary  layer  approximation.  The  ratio  ^  is  ssall  and  allows  the 

dissipation  terms  in  the  energy  equation  to  be  neglected  (14:78). 

All  other  terms  are  retained.  Finally,  assuming  the  variation  in  c 

p 

/  ^ 

due  to  turbulence  is  negligible,  the  terms  h  and  h  are  written  as  c  T 

P 

/ 

and  c  T  respectively.  The  system  of  equations  becomes: 

p 


Continuity 


a(rp)  a(rpu) 

-  +  -  =  0  (B.14) 

at  ax 


x-momentum 


a  __  a  _  _  a  _  _-r, 

— (rpu)  +  — [r(puu  +  p)]  +  — [r(-T  +  puv)]  =  0  (B.15) 

at  ax  ar 


Energy 


3  ^  5  _  3  _  _  -r,  _ap 

— (rpc  f)  -  — (rp)  +  — (rpuc  f)  +  — r(q  +  pc  vT)  -  ru —  =  0 

at  ”  at  ax  ar  ”  ax 


(B. 16) 


The  Navier-Stokes  equations  are  in  a  simplified  form  suitable  for 
deriving  the  damped-wave  equations. 
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The  continuity  equation  may  be  expanded  to  produce  the  following: 


dp  du  _dp 

—  +  p —  +  u —  =  0  (B.17) 

at  ax  ax 


Expressions  for  ^  and  c®"  be  obtained  from  the  equation  of  state: 


ap 

'i  ap  1  af 

—  =  p 

at 

p  at  T  at 

ap  _ 

'i  ap  1  af 

—  »  p 

ax 

p  ax  T  ax 

(B.18) 


(B. 19) 


Equations  (B.18)  and  (B.19)  may  be  written  in  a  more  compact  form  by 
recognizing  that  ^  =  d(Ln(pl)  : 


. 

'a(Ln[p])  1  af 

—  =  P 

- - - 

at 

at  T  at 

\  / 

ap 

'a(Ln(p])  1  af 

-  =  P 

3x 

ax  f  ax 

V  V 


(B.20) 


(B.21) 


Substitution  of  Eqs  (B.20)  and  (B.21)  into  the  the  continuity 
equation  produces  the  following  result: 


p 

■  ai  _aT  ■ 

au 

- 

-1 

— 

-  +  u — 

=  p —  +  p 

— Ln[p]  +  u — Ln[p] 

T 

_  at  ax 

ax 

[at 

ax  J 
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The  energy  equation,  Eq  (B.16),  can  be  written  as 


_ 'ST  _aT]  a(rp)  .dp  3  f  _  -7^,] 

rc  p  —  +  u —  -  -  -  ru —  +  — r  q  +  pc  vT  =0  ( B .  23 ) 

**  at  ax  at  ax  ar  ** 


or  more  compactly  to  produce 


_af 
a —  - 
ax 

4 


rc  p  -  +  u —  -p— r(Ln[p))  -  rup 


__a(Ln[p])  a 


-  +  — r  q  +  pc  vT  =0 

ax  ar  ^ 


(B.24) 


Equation  (B.22)  can  be  used  to  replace  the  first  term  of  Eq  (B.24): 


_ ’au  'a(Ln(p])  _a(Ln[p])y  _a  _  _ a(Ln[p)) 

rc  pT  —  +  -  +  u -  -  p — r(Ln(p])  -  rup - 


— r(q  +  pc  vT)  =  0  (B.25) 

ar  '■  P 


The  relations  p  =  pRT  ,R  =  c-c^,  used  to  simplify 

^  V 


Eq  (B.25): 


au  3(Ln[pl)  _  3(Ln(pl)  (y  -  1)  3  _  _  —r, 

ry —  +  r  -  +  ru  -  + - r(q  +  pc  vT)  =0  (B.26) 


p  dr 


r  p 


The  momentum  equation  Eq  (B.15)  is  rewritten  using  the  relations 
p  =  pRT  ,  and  c^  =  yRT 


au  c  a  _  3u  c  a  , 

r —  + - r(Ln[p])  +  ru —  +  ^ - [r(-T  +  puv)]  =  0  (B.27) 


at  y  ax 


3x  py  ar 
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The  desired  form  of  the  damped  wave  equations  is  obtained  by  isolating 
the  friction  and  heat  transfer  terms  on  the  right  hand  side  of  the 


equations: 


a(Ln[pl)  dxi  _3(Ln(pl)  (y  -  1)  1  3  _  -  -r, 

-  +  y —  +  u -  - - ^ - [r(q  +  pc  vT)]  (B.28) 

3t  3x  3x  p  r  3r 


3u  c^3(Ln[pl)  _  3u  c^  1  3  _  r, 

—  + - +  u  —  =  ^ - [r(T  -  puv)]  (B.29) 

at  y  ax  3x  py  r  ar 


Fig.  B. 1  Partitioning  of  Flow 


For  each  flow  interval,  averaged  temperature  and  averaged  velocity 
functions  are  defined  for  each  region  shown  in  Figure  (B.l): 


T  (x,t)  =  — rr-f  T(x,r,  t)2Trrdr 

n  AA  J 

n  r 

n-  I 
r 

ip” 

u  (x.t)  =  ■■■  ■  ■  u(x,  r,  t)27rrdr 

n  aA  J 

n  r 

n- 1 


(B.30) 


(B.31) 
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These  averaged  fxinctions  are  Incorporated  into  Eqs  (B.28)  and 


(B.29)  as  follows: 
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n 
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du 

— Zirrdr  + 

at 

1 
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- - 2irrdr  + 

AA  r  y  ax 
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’'r  ax 


1  Sn  a 
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n  n-1 


AA  r  py  ar 
n  n-1 


(B.32) 
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All  variables  are  time-averaged  except  u  and  v  .  Evaluation 
of  the  integrals  produces 


au  c^a(Ln[p])  a 

_"  + - +  _ 

at  y  ax  ax 


2nc^ 

- fr  (t 

ypAA  ^ 

n 


-  puv) 


I  9 

r  (t  -  puv) 
n-1  xr  r=r 

n- 1 


(B.33) 


2 

In  Eq  (B.33),  c  has  been  assumed  to  be  constant  over  the  line 
cross  section  and  is  taken  as  v'yRT.  Here  T  is  the  average  fluid 
temperature  taken  over  the  cross  section  «uid  y  =  1.4  .  Some  studies 
suggest  values  of  y  less  than  1.4  to  provide  better  agreement  with 
measured  results  since  these  smaller  values  of  y  account  for  heat 
transfer  to  the  pipe  walls  (19:11-12).  Heat  transfer  to  the  pipe  walls 
is  already  accounted  for  by  use  of  the  energy  equation  with  the  damped 
wave  equations.  As  a  result,  the  adiabatic  speed  of  sound  is  used  and 
y  is  set  at  1.4.  The  use  of  the  average  temperature  of  the  cross 
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section  to  determine  c  is  justified  from  the  assumption  of  uniform 
pressure  over  the  cross  section.  Also,  the  temp>erature  variation  over 
the  cross  section  away  from  the  wail  will  not  be  large. 


In  order  to  solve  for  u  in  Eq  (B.33],  the  assumption  that  (u  ) 

n  n 

2 

can  be  approximated  by  (u)  must  be  made.  For  a  uniform  velocity 

n 

profile  this  is  exactly  true  but  for  nonuniform  velocity  profiles,  this 
can  only  be  considered  as  an  approximation.  However,  the  velocity 
profiles  will  be  shown  to  be  relatively  flat  away  from  line  wall. 

Also,  the  flow  interval  near  the  wall,  where  the  largest  velocity 
gradients  exist,  contains  only  a  small  portion  of  the  flow  energy.  As 
a  result,  approximating  (u^)  by  (u)^  ,  should  not  introduce  a 

n  n 

significant  error  in  the  results  for  the  overall  flow.  The  equation 
for  averaged  velocity  function  may  now  be  written  as 


du  c^d(lnp)  3u  2itc^ 

n  n 

—  + - +  U  —  =  - 

n . 


at  r  ax 


ax  rpAA 


fr  (t 

L  n  xr 


-  puv) 


t  / 

r  (t  -  puv)  ] 

n-1  xr  r  =  r 


n*  1 


(B.34) 


In  a  similar  manner  the  equation  for  Ln(p)  is  obtained  in  terms  of 
averaged  variables: 


a(Ln[p] )  a(Ln(p] )  du 

-  +  u  -  +  7 — 

at  dx  dx 


27ra(r  -  1) 
- 

AA  p 


+  pvc  T) 

p  rsa 


(B.35) 


a  a 

where  ^  [  u27rrdr  and  AA  =  f  27rrdr. 

■'o  "^0 
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The  friction  term  on  the  right  hand  side  of  the  Eq  (B.34)  is 


written  as 

(B.36) 

The  heat  transfer  term  in  the  Eq  (B.35)  is  written  as: 


(B.37) 


The  variables  u  and  k  are  the  sum  of  the  molecular  and  turbulent 

T  T 

viscosity  and  heat  conduction  coefficients,  respectively.  Values  for  4^ 
are  found  using  a  suitable  turbulence  model  and  values  for  k^  are  based 
on  the  Reynolds  analogy. 

Neither  Eq  (B.34)  or  (B.35)  provides  a  meauns  for  finding  T^  ,  thus 
a  third  equation,  the  energy  equation,  is  required  to  find  T  .  The 

n 

energy  equation  Eq  (B.16)  is  rewritten  as 

aT  aT  c^  a  c^  a(Ln[p])  c^  a  r  aT 

r —  +  ru —  - r(Ln(p]) - ru  -  = - rk  — 

at  ax  yc  at  yc  ax  yc  p  ar  ^ar 

p  p  p  '• 

(B.38) 

Applying  the  definition  of  the  average  temperature  function  to  the  Eq 
(B.38)  produces 


aT  aT  c^  a  a(Ln{p] ) 

— "  +  u  — - (Lnlpl) - u  -  = 

at  "ax  yc  at  yc  "  ax 

p  p 


,2  2  . 

yc  p  (r  -  r  ) 

p  n  n-1 


r  k  (r  ) 

n  T  n 


-  r  k^(r  ) 

n-1  T  n-1 

r*r 

n 


(B.39) 
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a(Ln[p] ) 

The  terms  -  amd  u  are  obtained  from  solution  of  Eqs  (B.34)  and 

dx 

(B.35).  The  term  (u  T)  In  Eq  (B.39)  is  assumed  to  be  approximated  by 

n 

2 

u  T  with  the  same  reasoning  used  in  the  approximation  of  (u  )  in  Eq 

n  n  n 

(B.34).  This  term  is  representative  of  the  average  amount  of  enthalpy 
transported  by  convection  in  the  interval  I 
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APPENDIX  C 


Velocity  Profile  Approximation  by  Piecewise  Continuous  Polynomials 


Piecewise  continuous  polynomials  are  used  in  Chapter  II  to  reduce 
the  number  of  grid  points  required  in  the  radial  direction  to  obtain  an 
accurate  nximerical  solution  for  transient  flow  in  a  fluid  transmission 
line.  The  polynomials  used  to  approximate  the  yelocity  profile  in  each 
of  the  radial  interyals  dlyiding  the  flow  are  shown  in  Eq  (C.l) 

(70:591): 


Pj(x,r,t) 

2 

ssa  +  br  +  cr 
11  1 

r  €  I 

1 

p  (x,r,t) 

n 

=  a  +  b  fr-r  ) 

n  n  '•  n- 1' 

+  c  fr-r 

n l  n- I' 

+ 

d  (r-r  ) ^ 

n'  n-1’' 

r  €  I 

n 

(C.l) 

For  n  =  1  • • •  M  ,  this  system  of  equations  contains  4M-1  unknown 
coefficients.  These  coefficients  must  be  evaluated  at  each  time  step 
and  axial  station. 

Application  of  the  definition  of  the  averaged  velocity  function  to 
Eq  (C.l)  produces  M  equations: 


u^ (x, t) 
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(C.3a) 
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r  5 

^  n-1  n 

+ 

o 

9 

3 

1  ^ 

2 

2 

r  5 

n-1  n 

3  r  +  r 

r  +  r 

2  r  +  r 

3 

r  +  r 

n  n-1 

n  n-1 

n  n-1 

n  n-1 

3  1 


1  r  5 

1  n-l  n 


5  r  +  r 

n  n-1 


2  r  +  r 

n  n-1 


n  =  1 


M 


(C.3b) 


The  piecewise  defined  polynomials  as  well  as  their  first  and  second 
derivatives  are  required  to  be  continuous  at  the  boundary  between  each 
flow  Interval.  These  conditions  produce  an  additional  3M-3  equations: 


(C.4) 
(C.5) 
(C.6) 

The  assumption  of  no  slip  at  the  wall  produces 


a  =a  +b5+c  5^+  d  5^ 

n>l  n  n  n  n  n  n  n 

b  =b  +2c3+3d5^ 

n^l  n  n  n  n  n 

2c  =  2c  +  6d  5 

n*l  n  n  n 


0=a  +b(r-r  )+c(r-r  )^+  d  (r  -  r  (C.7) 

M  MM  M-1  M  M  M-1  M  M  M-1 


euid  the  assumption  of  axisymmetric  parallel  flow  _^=  0  is  applied 
at  the  centerline: 


b^  =  0  (C.8) 

producing  two  more  conditions.  Equations  (C.3)  through  (C.8)  provide  a 
system  of  4M-1  independent  equations  to  solve  for  the  4M-1  coefficients 
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In  Eq  (C.l). 

The  above  equations  need  not  be  solved  as  they  appear,  but  may  be 
used  to  provide  a  set  of  linear  dependence  relationships  (70:592). 

Note  that  d^  =  0  since  the  equation  in  the  center  interval  is 
quadratic.  With  d^  *  0  ,  Eq  (C.6)  produces  the  relation  c^  =  c^ 

when  n  »  1  Equation  (C.6)  is  solved  for  d  to  give 

n 

-(c  -  c  )  n  =  2,3.---,M-l  (C.9) 

n 

By  substituting  Eq  (C.8)  into  Eq  (C.5)  an  expression  for  b  is 
produced: 


d  = 

n 


b  =*b+5(c+c  )  n=l,2,'-’,M  (C.IO) 

n>l  n  n  n  n^l 


For  specific  values  of  n,  Eq  (C.IO)  can  be  written  as 


2 

b 

3 
b 

4 


=  0  +  5  (c  +  c  ) 

11  2 

=b  +5(c+c)  =5(c+c)  +5(c2+c3) 

2  2  2  3  1  1  2  2 

=b  +5(c+c)=5(c+c)+5(c+c)+5  (c3  +  c4) 

3  334  112  223  3 


(C. 11) 


A  recursion  relationship  for  b  in  terms  of  c  can  be  developed  from 
Eq  (C.ll)  by  inspection: 
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n 


n*X 


(C. 12a) 


b 

n-*’  1 


C 

k  k 


C 

k 


k=l  k=2 


or  equivalently 


n- 1  n 

b  =  y  5  c  +  y  5  c  n  =  2,  3. • • • ,M  (C. 12b) 

n  L  k  k  L  k-1  k 

k=l  k=2 


Equation  (C.4)  may  also  be  written  for  specific  values  of  n.  In 
addition,  this  result  may  be  written  in  terms  of  c^  only,  if  Eqs  (C.9) 
and  (C.ll)  are  used  to  eliminate  the  coefficients  b  and  d  : 

n  n 


a  =  a+  b5+c  d  5^  *  a  +  c  5^ 
2  1  11  11  11  1  11 


(b^=0) 


(d  =0) 
1 


=  a  +  b  3-+  c  d  3^ 
2  2  2  2  2  2  2 


=  a  +  c  3“+  3  [3  (c  +  c  )]+  3‘‘c  3^r— (c  -  c  )! 

1  1  1  2  1  1  2  2  2  2l-  3  2-1 

Jo 

2 

a  =  a  +  b3+c  3^+  d  3^ 

4  3  33  33  33 

=  a  +  c  3^+  3  [3  (c  +  c  )]  +  5^c  +  -  5^(c  -  c  )  + 

1  11  211  2  22^23  2 

3  [5  (c  +  c  )  +  3  (c  +  c  )]  +  3^c  +  3^r~(c  -  c  )] 

3  1  1  2  2  2  3  3  3  3*-  4  3 

Jo 

3 

=  a-f3(3c)  +  5  (3  c)  +  5  (3  c)  +  3  (3  c) 

1  111  211  211  222 

5(3c)  +  3  (3  c)  ■'•o(3c)  +  5  (3  c)  + 

311  312  322  323 


3^  3^ 

3  (3  c  )  +  — (c  -  c  )  -t-  — (c  -  c  1 

3  3  3  332  3^3 


(C. 13) 


By  inspection  of  Eq  (C.13),  the  pattern  formed  by  the  3  and  c 

n  n 
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coefficients  produces  a  general  expression  for  a  : 


n-l  J  n-l  J  n-l 


Y  d  ] 

r  5  c  +  y  5  ] 

(5  c  +  —  y3^(c  -c) 

(C.14) 

L  J  L 

L  k  k  L  J  ( 

L  k-1  k  L  J  j^i  J 

J=1  k=l  J=2  k=2  J=2 


Note  that  a  is  in  terms  of  S  and  c  except  for  the  term  a  . 

n  n  n  1 

Substitution  of  Eq  (C.7)'for  a^  on  the  right-hand  side  of  Eq 
(C.14)  evaluated  at  n  =  M  results  in  the  following: 


M-l  J  M-l  J  M-l 


a  +  ] 

r  3  ] 

1  3  c  +  ^ 

'  3  ] 

^  It 

(3  c  +  -  ) 

'  3^(c  -  c  ) 

+ 

1  L 

-  J  1 

^  k  k  1 

-  J  !■ 

-  k-1  k  L 

-  J  J 

J=1  k=l  J=2  k=2  J=2 


M-l 


I* 


c  + 

k  k 


N 


C 

k-1  k 


k=«l 


ks2 


+  5^c  +  5^d  =  0 

MM  MM 


(C.15) 


Equation  (C.15)  is  solved  for  a^  to  produce 


a  =  - 
1 


M  J  M  J  M-l 

1  r  .2, 


Yd  c  +  iy  5^(c  -  c  )  +  5^d 

L  J  L  k  k  L  J  L  k-1  k  2  L  J  J*l  J  Ml* 

j  =  l  k  =  l  J  =  2  k=2  '^J=2 


(C. 16) 


Substitution  of  Eq  (C.16)  into  Eq  (C.14)  produces 


a 

n 


Y  5  Y  5  c  - 
L  J  L  k  k 


J=1  k=l 


M-l  J 

Yd  Ydc- 

L  j  L  k  k 


J=1  k=l 


M-l 

5  Y  d  c  -d^c 

H  L  k  k  H  H 
k  =  l 


+ 
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n-1  J  K-1  J  M 


7  5  c  - 

'a  c -  5  1 

)  S  C 

+ 

L  J 

L  k-l  k 

L  J  i 

L  k-l  k  M  ^ 

Li  k-l  k 

J»2  k=2  J=2  k=2  k=2 


C 


3^d  (C.17) 

M  M 


The  terms  In  the  brackets  may  be  combined  to  produce 


H  j  M  J  M-l 


a  =  “  ] 

r  5  ] 

r  5  c  - ; 

r  5  ] 

: s  c-i) 

r  5^(c  -  c  )  -  5^d 

n  L 

-  J  i 

k  k  i 

-  J 1 

-  k-l  k  _  L 

-  J  J+l  J  MM 

J=n  k=l  J=n  k=2  J=n 


(C. 18) 


The  recursion  relationships  for  all  the  unknown  coefficients  of  the 
piecewise  defined  polynomials  of  are  summarized  below: 


M  J  M  J  M-l 


a  =  “  ] 

)  5  ] 

r  5  c  -  ] 

r  5  ] 

r  a  c  -  ‘ ; 

'  5^(c  -  c  )  -  5^d 

n  i 

-  J  !■ 

k  k  / 

-  J  t. 

-  k-l  k  L 

-  J  J*1  J  MM 

j=n  k=l  J=n  k=2  J=n 


n  =  1,  2,  •  •  • ,  M 


b  =  0 
1 

n-1  n 

b  =  750+^5  c 

n  L  kk  L  k-lk 

k=l  k=2 


d 

1 


=  0 


d 

n 


— ^(c  -  c  ) 

n>l  n 


n  =  2,  3, • •  •  ,M 


n  =  2,  3. • • • ,  M-l 


(C. 18) 
(C.8) 

(C. 12b) 

(C. 19) 
(C.9) 


All  polynomial  coefficients  are  now  expressed  In  terms  of 
c  ,c  ,...c  auid  d  ,  which  must  be  determined  from  values  of  the 

2  3  H  M 

averaged  velocities.  This  is  accomplished  using  Eqs  (C.2)  and  (C.3) 
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which  describe  the  averaged  velocity  in  each  interval: 


where 


a  +  -5b  +  -5c  =  u 

1  3  '  '  2  '  '  ' 

a  +  m  b  +  m  c^+  m  d,  = 

2b2c2d2  2 

2  2  2 


a  +  m  b  +  m  c  +  m  d  =  u 

n  bn  c  n  dn  n  n 
n  n 


m  -  -  8 
^  3  ^ 

m  *  i  5^ 

S  2  ^ 


m  =0 
d 

2  n  n 

m  - - +  r  - 

b  _  n-l 

n  3r  +  r  r  +  r 

n  n-1  n  n-1 

1  3^  8^ 

I  n  d  n 

ra  - - +  -  r  - 

'^n  2  r  +  r  3  r  +  r 

n  n-1  n  n-1 

2  3^  1 

2  n  1  n 

m  = - +  -  r  - 

d  _  _  n-1 

n  5r  +  r  2  r  +  r 

n  n-1  n  n-1 


(C.20) 


(C.21) 

Notice  that  coefficients  m  ,  m  ,  and  m  in  Eq  (C.21)  depend 

be  d 

n  n  n 

only  on  the  known  values  of  r  through  r  which  are  the  locations  of 

1  n 

the  flow  intervals.  These  values  of  r  through  r  are  selected  at  the 

1  n 


beginning  and  never  change  over  the  course  of  the  solution.  Also,  the 
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values  of  are  considered  known  since  they  are  provided  by  the 
solution  of  the  damped  wave  equations  or  as  part  of  the  set  of  initial 
conditions. 

The  coefficients  a^,  b^,  and  d^  in  Eq  (C.20)  can  be  expressed  in 
terms  of  *  circular  line  with  four  intervals, 

Eqs  (C.2)  and  (C.3)  may  be  written  as 


d^  =  u 

4  4  1 


m 


c  [23  m  +m  - — -~ZS  (6  *S  *6  )-3  (5  *6  *6  + 

^  '  ‘*2  '2  33^  1  2  3  4  223432 


m 


c_t— -5_(3  +3  )-3  (3  +3  )-i(3^-3^)]+c  [-3  (5  +5  )-V]  + 

3  334  234  23  4  434  3 

2 


d  (-3  )  =  u 

4  4  2 


c^[m^  (23  +3  )-23  (5  +5  )-3  (3+5  )]+ 

2b^  12  134  234 


m 

d 

C^[m^  5  +m - --6  (3  +5  )-5  (5  +5  )+^^]  + 

3b2c  33d  P  ld  T 

3  3  35  3 

3 


m 


c  [~-5"-5  3  -i3‘]+d  (-5^)  =  u 

4  35  44333  4  4  3 

3 


(C.22a) 


(C.22b) 


(C.22c) 
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(C.22d) 


c  Im  5  +in  -(3  +5  5  )]+d  [m 

4  b3c  443  4d 

4  4 


U 

4 


In  the  above  equations,  all  expressions  in  the  square  brackets  are  in 
terms  of  3  through  3  ,  or  equivalently,  r  through  r  .  The  value  of 

14  14 

these  bracketed  terms  are  known  and  behave  as  constaints  once  the 

locations  of  the  interval  boundaries  have  been  selected.  Let  1  ,  m  , 

1  1 

p^,  and  q^  represent  the  bracketed  terms  in  Eq  (C.22)  which  multiply 
c  ,  c  ,  c  ,  and  d  respectively.  Equations  (C.22a-d)  are  now  written 

2  3  4  4 

in  matrix  form  as 


m  p  q 
1  1  1 


m  p  q 
2  *^2  ^2 


m  p  q 
3  *^3  ^3 


m  p  q 

4  ^4  ^4 


c 

U 

2 

1 

C 

u 

3 

2 

C 

U 

4 

3 

d 

U 

J 

4 

4 

(C.23; 


The  unknowns,  c  ,  c  ,  c  ,  and  d  ,  are  obtained  by  solving  Eq  (C.23): 

2  3  4  4 


-1 

C 

r  1 

m 

p 

q 
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2 

1 

1 

1 

1 

1 

C 

1 

ni 
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q 

u 
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*^2 

^2 

2 

C 

1 

m 

p 

q 
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4 

3 

3 

^^3 

^3 

3 

d 

1 

m 

P 

q 

U 

4 

4 

4 

*^4 

^4 

4 

(C.24) 


This  velocity  profile  approximation  may  be  extended  to  any  number  of 
intervals.  The  method  for  determining  the  unknown  coefficients  in  the 
temperature  profile  polynomials  is  exactly  the  same  as  shown  above. 


229 


The  profile  of  T  -  is  used  Instead  of  T  to  produce  a  value  of  zero 
at  the  wall. 
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APPENDIX  D 


Unsteady  Turbulent  Friction  Model 


From  an  analogy  with  laminar  flow,  the  eddy  viscosity  formulation 
of  Bousslnesq  has  the  form  (61:180): 


au 

T  =  -  puv  =  ^ 

Turb  ^  Turb  dy 


Prandtl  suggested  the  mixing  length  formulation: 


(D.  1) 


T  =  -puv  =  pt 

Turb  ^  m 


au 

ay 


au 

ay 


(D.2) 


where  t  Is  some  effective  Interaction  dlstauice  between  eddies.  A 

■ 

relationship  for  i  as  a  function  of  the  dependent  variables  Is 

■ 

required  to  calculate  t  The  eddy  viscosity  and  the  mixing  length 

Turb 

formulations  are  essentially  equivalent  since  inspection  of  Eqs  (D.l) 
and  (D.2)  shows 


^Turb 


ay 


(D.3) 


The  steady  state  model  for  apparent  viscosity  may  be  used  as  a 

starting  point  in  developing  a  model  capable  of  describing  transient 

effects.  Very  near  the  wall,  viscous  shear  dominates  producing  an 

approximately  linear  velocity  profile:  t  =  p  -  for  0  s  y*  s  5  and 

w  y 
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y*  ,  where  = 


yu. 


u  . 

u  *  -  ,  and  u.= 


.  Further  f rom 


V'  u.  ■  •  V  p 

the  wall  both  viscous  and  turbulent  shear  become  Importamt  but  the  flow 
is  still  dominated  by  the  pipe  wall.  The  velocity  profile  in  this 
logarithmic  overlap  region  found  at  50  s  y*  s  500  is 


u*  =  ^  In(y^)  +  C  (D.4) 

The  region  described  by  5a  y*a  70  may  be  considered  as  a  transition 
region  over  which  an  expression  must  connect  the  velocity  profiles  from 
the  viscous  sublayer  and  the  inner  layer.  The  linear  velocity  profile 
of  the  viscous  sublayer  must  blend  smoothly  into  the  logarithmic 
velocity  profile  of  the  inner  part  of  the  inner  layer  (61:182). 

Van  Driest  developed  a  means  of  describing  the  shear  stress  in 
this  transition  region  (61:183).  The  shear  stress  cam  be  expressed  in 
the  transition  region  as: 


T  +  T 


Turb 


du 


puv  =  T 


(D.5) 


where  t  is  the  laminar  shear  stress.  In  the  logarithmic  part  of  this 

tramsition  region  the  turbulent  shear  stress  dominates  so  that 

T  »  T  and  T  a  T  Applying  the  eddy  viscosity  formulation 

Turb  W  Turb 

to  T  ,  Ed  (D.l)  can  be  written  as 

Turb 


du 

^Turb  dy 


T 

H 


(D.6) 
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Additionally,  the  derivative  of  the  log  law  in  this  region  produces 


^  _  _1 

dy  ~  icy 


(D.7) 


Substituting  Eq  (D.7)  into  Eq  (D.6)  produces  an  expression  for  the  eddy 
viscosity: 


^«rb  =  «  p  u.  y  »  P 

Vaui  Driest  represented  the  damping  of  turbulent  eddies  due  to  the 

presence  of  a  rigid  wall  by  an  oscillating  wall  in  a  fluid  at  rest. 

This  problem  had  previously  been  solved  by  Stokes  (61:93),  Van  Driest 

applied  the  factor  1  -  exp(-  ^),  based  on  Stokes’  solution,  to  obtain 

A 

the  damping  effect  of  the  wall.  The  result  of  Eqs  (D.6)  and  (D.8)  is 
used  with  Eq  (D.5)  to  produce 


T  +  T 


Turb 


du  ^  2  af 

1  f  “y  1 

du 

=  =  p  .  p  x  y  ^ 

l-exp^  -X  J  J 

dy 

au 

dy 


^  (D.9) 


Van  Driest  found  that  A  =  —  and  ic  =  0.40  accurately  represented 

data  in  the  wall  region. 

The  Van  Driest  mixing  length  model  for  the  entire  inner  region  is 


t  =  K  1 


exp 


26v 


y 


(D. 10) 
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An  expression  for  the  eddy  viscosity  is  found  by  comparison  of  Eq 
(D.IO)  with  Eqs  (D.l)  and  (D.2): 


2 

■ 

'  -yu. 

1^  2 

au 

»  p  1C 

1  -  exp 

26v 

C  J 

y 

ay 

This  eddy  viscosity  model  is  valid  only  for  fully-developed, 
steady,  turbulent  flow.  The  effect  of  the  unsteadiness  resulting  from 
transient  flow  disturbances  produced  by  traveling  pressure  waves  must 
be  included  in  this  model  to  make  it  suitable  for  the  purposes  of  this 
study.  Cebecl  accounts  for  unsteady  effects  by  using  instead  of 
in  the  calculation  of  u^  (12:2153).  The  variable  represents  the 
shear  stress  evaluated  at  the  edge  of  the  viscous  sublayer  rather  them 
at  the  wall.  This  is  Justified  on  the  physical  grounds  that  the  random 
velocity  fluctuations,  damped  according  to  the  mechanism  of  the  Stokes 
layer,  do  not  extend  all  the  way  to  the  wall  but  disappear  at  the 
viscous  sublayer.  Cebeci  approximated  the  momentiim  equation  in  the 
immediate  vicinity  of  the  wall  of  a  flat  plate  to  include  only  the 
effect  of  the  pressure  gradient  (12:2153): 


dx  ap 

s  ^ 

ay  ax 


(D. 12) 


Telionis  recommends  the  unsteady  velocity  term  be 
approximation  of  the  momentum  equation  by  modifying  Eq 


included  in  the 
(D.12)  (66:233): 


au 

ap 

dx 

s 

—  = 

—  —  ■¥ 

at 

ax 

ay 
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(D. 13) 


T  Is  found  by  integrating  from  the  wall  to  to  the  edge  of  the  viscous 
sublayer  denoted  by  y*  .  Cebecl  chose  y^  =  11.8  as  the  boundary 
between  the  viscous  sublayer  and  the  log  law  layer  based  on  the 
intersection  of  the  log-law  equation  and  the  equation  of  the  viscous 
sublayer  for  steady  flow.  The  expressions  \i*  ^  -  In(y^)  +  C  and 

K 

u*^  »  y*  are  equal  at  y*  *  y*  The  intersection  of  the  two 
equations  can  easily  be  seen  In  the  universal  wall  law  plot  for  a 
turbulent  boundary  layer  on  a  smooth,  solid  surface  (73:476). 

Transforming  Eq  (D.13)  to  cylindrical  coordinates  for  pipe  flow 
and  Integrating  from  the  edge  of  the  viscous  sublayer  to  the  wall 
produces 


a(rT) 

dr 

r 

a 


dr 


au 

at 


r  dr  + 


dx 


r  dr 


Note  that  the  coordinate  y  in  Eq  (D.13),  measured  out  from  the  wall, 
has  been  replaced  by  the  coordinate  r  whose  origin  is  now  at  the  center 
of  the  fluid  line.  Assuming  that  p(p,T^)  =  near  the  wall, 
integration  of  the  above  equation  produces: 


ax  - 


r  T 

a  a 


=  P, 


au  ,  , 

^  r  dr  +  (a 


-  r 


)  ^ 
'  dx 


or  solving  for  t  : 

S 


235 


(D. 14) 


T 

s 


ax  +  (r  -  a)  ^  +  p 
w  ■  dx 


3vi  . 
at  " 


The  term,  p 


au 

—  r  dr 

at 


in  Eq  (D. 14)  must  be  evaluated. 


The  inner  boundary,  of  the  region  containing  the  viscous 

sublayer  is  chosen  to  be  less  than  r  ,  where  r  =  a  -  y  and 

vy^  ■  s  s 

=  — -  The  velocity  profile  in  this  region  is  approximated  by  the 

piecewise  continuous  polynomials  as 


u(r)  =  a  +  b  (r  -  r  )  +  c  (r  -  r  )^+  d  (r  -  r  )^ 

H  M  M-1  M  M-1  M  M-1 

(D.15) 

Equation  (D.15)  is  the  cubic  velocity  profile  approximation  used  in 
Chapter  II.  For  r^  contained  in  the  interval  adjacent  to  the  wall, 
integration  of  this  velocity  profile  function  from  r  to  the  wall 

S 

produces 


r  Z 


-[(r  -r  )^-(r  -r  )^] 
3L  M  M-1  s  M-1  J 


^  r  il 
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Equation  (D. 14)  cam  now  be  evaluated  at  each  time  step  as  the  equations 
developed  In  Chapter  II  and  III  are  solved  numerically.  Once  Is 
found,  the  unsteady  value  for  the  Vam  Driest  damping  factor  A  Is  easily 
calculated  at  each  time  step: 


A 


261/ 

u 

S 


(D.17) 


where  ®  ^  expression  for  the  mixing  length  for  pipe 

flow  cam  be  written  from  Eq  (D.IO): 


I  =  0.41(a 


r)a 


-  exp 


(D. 18) 


This  expression  for  the  mixing  length  is  an  extension  of  the  model 
developed  by  Cebeci  in  which  Van  Driest’ s  damping  factor  is  modified  to 
include  the  effect  of  the  time  varying  velocity  in  addition  to  the  time 
varying  pressure  gradient  (65:1000). 

The  Clauser  eddy  viscosity  model  is  assumed  to  hold  in  the  outer 
layer  for  flat  plate  flow: 


c  =  p  K  U  5 

o  2  c 


where  =  0.0168  and  5  is  the  displacement  thickness.  The  eddy 
viscosity  is  assiamed  to  be  constant  across  the  outer  layer.  White 
points  out  that  the  outer  layer  is  commonly  neglected  in  pipe  flow 
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since  it  has  little  effect  on  tho  analysis  and  the  lack  of 
intermit tence  in  duct  flow  produces  greatly  increased  turbulence  at  the 
centerline  as  compared  to  flow  over  a  flat  plate  (73:473).  As  a 
result,  Eq  (D.ll)  will  be  used  to  calculate  eddy  viscosity  from  the 
wall  to  the  centerline  of  the  transmission  line  with  no  special 
expression  governing  the  outer  layer. 


238 


APPENDIX  E 


Development  of  Modal  Superposition 


Separation  of  Variables 

To  find  the  modal  representation  of  Eq  (90)  f  is  expressed  as  the 
product  of  mode  shapes  G^{x)  times  modal  general  coordinates  ^^(t) 
(31:429).  Thus  f  is  written  as 


f  =  EG.(x)  ?  (t) 
1=0 


(E.  1) 


This  technique  is  commonly  applied  to  linear  partial  differential 
equations  and  is  known  as  separation  of  variables,  the  product  method, 
or  Fourier’s  method  (53:50).  Homogeneous  boundary  conditions  are 
required  for  this  method  to  work.  Substitution  of  Eq  (E.l)  into  Eq 
(90)  produces 


p 

PF 
+  - 

2  2 

pc  a  G 

at^ 

at 

AG  ax^ 
1 

The  forcing  function  on  the  right-hand  side  of  Eq  (88)  has  been  set  to 

zero  so  that  the  modal  frequencies  can  be  found  (32:82).  In  the 

3f* 

absence  of  a  forcing  function,  -5—  must  vanish  at  each  end  of  the  line 

ax 

as  a  result  of  Eq  (87).  Since  and  —  are  zero,  separation 

ot J  x=o  atj x=L 

of  variables  requires  the  following  to  be  true: 
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3\ 


at 


+ 


at 


G  ax^ 


2 


(i) 


(E.3) 


Equation  (E.3)  provides  an  expression  for  the  mode  shapes: 


3^G 

- -  +  G  w  =  0 

ax"  ‘  ‘ 


(E.4) 


The  boundary  conditions  for  Eq  (E.4)  are 

(E.5a) 
(E.Sb) 

A  general  solution  is  assumed  for  G^(x): 


g(O.t)  =  0 


ga.t)  .0 


G^(x)  =  B  cos(^)  +  D  sin(^) 


(E.6) 


The  boundary  conditions  of  Eq  (E.5)  require  that  D  =  0.  Assuming 
B  =  1  the  substitution  of  Eq  (E.6)  into  Eq  (E.2)  auid  Eq  (E.2)  into  Eq 
(E.l)  produces 


P  ,  5"?  pF  .  a^  pc"i"7r"  . 

,i7rx.  1  /iJrx-i  i  <- 

-  COS(-; )  -  +  COS(— ; )  +  -  COS(-;^ — )  t 

A  at"  A  ^  at  al"  ' 


e  +  e 
1  2 


(E.7) 


The  orthogonality  of  the  mode  shapes  may  be  used  by,  first,  multiplying 
both  sides  of  Eq  (E.7)  by  cos(^^^)  and  then  integrating  over  x  from  0 
to  L.  The  dependence  on  x  is  removed  from  Eq  (E.7): 
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(E.8) 


pL  PLF  aC.  pc^i^ir^  .  . 

- i  + - ‘+  -  ?  =  e  +  e  (-1)‘ 

2A  at  2A  at  2AL 


Integration  over  time  produces 


pL  af  pLF 

- +  -  C  =  ®2 

2k  at  2k 


(E.9) 


pL  ac,  pLF  pc^iV^t 

- L  +  -  ^  +  - 


2A  at  2k 


1 


2AL 


I  ?  dt  =  e  +  e  (-1)*  (E.IO) 

J  _  1  i  2 


Karnopp  points  out  that  Eq  (E.9)  represents  the  line  as  a  single, 

lumped  Inertia  while  Eq  (E.IO)  represents  a  collection  of  modal 

oscillators  (31:439). 

Residual  Compliance 

The  residual  compliance  is  calculated  from  Eq  (E.IO)  at  steady 
state  with  the  mass  and  resistance  elements  suppressed  (32:84).  First, 

define  two  new  variables;  q^  =  J  ^^dt  amd  Pj  =  J  where  q^  and 

p^  represent  the  generalized  displacement  and  the  generalized  momentum 
of  the  ith  mode,  respectively.  Equations  (E.9)  and  (E.IO)  may  be 
rewritten  in  terms  of  these  new  variables; 

pL  a^q  pLF  aq 

- r  + - =  e  +  e  (E.lla) 

2k  at  2k  at  12 
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pL  S^q  pLF  dq  pc^l^n^ 

- -  + - -  *  -  q  =  e  +  e  (-1)  (E.  lib) 

2A  at^  2A  at  2AL  ‘  ^  ^ 


(E.llb) 

dominates  the  transient  response: 


pc 


AL  r  ^ 

-  e  +  e  (-1 ) ' 

,2  2  I  1  2  I 

i  TT  I-  J 


2AL 
~2 


(E. 12) 


This  expression  may  be  thought  of  in  terms  of  the  displacement  of  a 
spring  due  to  the  application  of  a  force  displacement  = 
force/stiffness.  Using  Eqs  (E. 1),  (E.6)  and  the  definition  of  q^  where 
q^  =  X^Cjdt  ,  Eq  (h. i2)  produces  the  following: 


-t  08  OB  2AL  r  1 

J  f,,  dt  >  IG,(0)  q^(t)  .  I  (1)‘ 

1=0  1=0  pc  i  71  >-  ■* 


pt  00  00  2AL  r  -J 

J  =  I  =,'■->  <1,“'  =  j:  I 

1=0  1=0  pc  i  71  >-  ■* 


(E. 13) 


Equation  (E.13)  may  be  written  in  the  following  form: 
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(E. 14] 
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Inspection  of  Eq  (£.14)  produces  the  residual  compliance  matrix.  Since 
this  residual  compliance  will  be  used  to  represent  modes  with 


frequencies  greater  than  ,  let  1  in  Eq  (E. 14)  become  n+J  where  n  is 
the  highest  mode  retained.  The  residual  compliance  is  now  written  as 


•  €D 

I 

J=i 

1 

o 

E 

j=i 

"c  c 

11  12 

2AL 

(!»♦  J) 

(-1)‘ 
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2  2 
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E 

j=i 
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21  22 

{ n+  J  ) 

Complementary  Equations 

The  same  elect'-ical  analogies  used  in  developing  Eq  (106),  (107), 
(111),  and  (112)  are  used  to  develop  a  complementary  system  of 
equations.  Development  of  this  system  of  equations  is  summarized 
below: 


Inputs  at  Each  End 


Fig.  E.2.  Bond  Graph  for  Equations  (E.16)  and  (E.17)  with  Mixed 
Pressure  and  Flow  Inputs  Modal  Equations  for  Both  Figures  E. 1  amd  E-2 


p  =  K  q 


q  =  -  M’^p  -  IR"\  q 

e 


(E. 16) 


(E. 17) 


Output  Relations  for  Figure  E.l; 


e  =  IH  p  +  11  f 

s  eq  e 


(E.18) 


Output  Relations  for  Figure  E.2: 


f  =-7^rre  dt-LHpl-^f 
si  1  I  J  el  1  1  02 

22  ^  J  22 


(E. 19) 


e  = 

s2 


1  - 
22 


1  1 
12  ' 


1  02  I  2  1  1 

22  •'  22 


■  ^  21 
p  +  —  e 

1  ol 
22 


(E.20) 
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The  variables  i  are  elements  of  the  residual  inertance  0  which  is 

l  i  eq 

obtained  in  the  same  manner  as  the  residual  compliauice  C 

eq 


0 

eq 


2pL 


OB 

z 

J=1 

1 

OB 

1 

( n*  J ) 

( n+  J  ) 

OB 

E 

j=i 

OB 

1 

J  =  i 

1  2 

(E 

(n+  J ) 

(n*J) 

The  elements  of  the  R,  K,  and  W  matrices  are  calculated  as  follows: 


Table  E. 1 

Inertia,  Capacitance,  and  Resistance  Parameters 


I 

C 

R 

1 

1 

i 

AL 

LAF 

1=0 

QB 

2 

2 

pc 

pc 

1=1, 2, 

2pL 

AL 

LAF 

2. 

Ink 

2pc^ 

2pc^ 

As  shown  in  Chapter  V,  it  is  possible  to  more  accurately  determine 
the  M,  K,  and  R  matrices  as  well  as  a  G  matrix  to  account  for  heat 
transfer  effects.  Equations  (E.19a,b)  are  written  in  terms  of 
frequency  dependent  parameters: 

a^q  aq 

L  I  +  L  R  —  =  0  (E.22a) 

at  at 
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L  d\  L  aq  L 
-  I  - i  +  -  R  — 1  +  - 

2  at^  2  at  2 


1  R  G  1 


(E.22b) 


The  complementary  formulation  for  Eqs  (E.17)  and  (E.IS)  are  written  as 


a^q  R  C  aq 

L  C  - -  +  L  (G  +  — )  —  =  0 

“  at^  °  I  at 
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(E.23a) 
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R  C 
1  1 . 


I 
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1  R  G  'i 


I 


I 


=  0 


(E.23b) 


Tables  E.2  and  E.3  replace  Tables  3  and  E. 1  respectively  for  Eqs 
(E.22a,b)  and  (E.23a,b). 


Modal  Frequencies 

The  modal  frequencies  for  the  systems  shown  in  Figures  3,  4,  E. 1, 
and  E.2  may  be  approximated  from  the  two  systems  of  modal  equations  and 
output  relations  developed  above.  The  trauismission  line  configuration 
shown  in  Figure  12  is  used  as  an  example.  Neglecting  the  damping 
terms,  Eqs  (106)  and  (107)  are  combined  and  written  as 

M  -  +  Kq-iH[e  .e  ]=0  (E.24) 

, . 2  el  s  2 
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Table  E.2 


With 


Table  E.3 


e  =0  and  f  =0  ,  Ea  (5.30)  produces 

el  e2 


e 

s2 


q 


Substitution  of  Eq  (E.25)  into  Eq  (E.24)  with  e  =0  produces 


(E.26) 


-  +  q  -  H^[0,  -  -L  L  H  qll  =0 

dt^  ^22  ^  } 

Assuming  a  solution  of  the  form  =  sin(«^  +  0)  for  Eq  (E.26) 
produces 


[k  +  L 
I  2  C 

V.  OO 


L 

2 


(E.27) 


where  I  is  the  identity  matrix  in  this  case.  The  modal  frequencies  may 
be  determined  by  setting  the  determinant  of  the  bracketed  term  of  Eq 
(E.27)  to  zero  and  solving  for  .  Similar  expressions  may  be 
developed  for  other  transmission  line  configurations. 
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